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PREFACE 


Precision  radar  and  optical  systems  require  very  accurate  tropospheric  re¬ 
fraction  corrections  to  account  for  time  delay  and/or  ray  path  bending  as  a  func¬ 
tion  of  the  atmospheric  state.  This  technical  note  describes  a  state-of-the-art 
raytrace  model  developed  and  operated  at  the  US  Air  Force  Environmental  Technical 
Applications  Center  (USAFETAC)  that  uses  a  single  vertical  meteorological  profile 
and  the  general  theory  of  geometric  optics  to  compute  atmospheric  refractive  ef¬ 
fects.  The  most  common  application  of  this  model  is  point-to-point  ray  traces. 
An  important  feature  of  this  model  is  its  ability  to  have  the  source  at  a  higher 
altitude  than  the  target.  In  other  words,  zenith  ray  angles  exceeding  90°  are 
possible.  The  output  format  of  this  model  features  a  table  with  values  of  errors 
due  to  refraction.  Actual  ray  plotting  is  not  provided. 

Most  readers  of  this  technical  report  will  not  want  or  need  to  read  the 
entire  report.  The  objective  will  probably  be  to  seek  out  specific  answers  to 
RAYTRA  questions.  The  Chapter  titles  will  help  the  reader  find  what  he  is  after. 
In  general,  the  person  who  wants  to  use  RAYTRA,  or  thinks  he  does,  should  read 
Chapters  2,  7,  and  8,  and  Section  4.3.  If  more  detail  or  background  is  required, 
refer  to  the  remaining  chapters. 

The  main  raytracing  model,  called  RAYTRA,  was  designed  primarily  by  Maj  John 
Mill.  The  programming  of  the  raytracing  package  (three  separate  programs,  one  of 
which  is  RAYTRA)  was  accomplished  primarily  by  Lt  Jack  D.  Brown,  Jr.,  with  tech¬ 
nical  help  from  Maj  Mill  and  Capt  Ted  Linn.  Capt  Michael  Abel  oversaw  the  pack¬ 
age  completion.  This  technical  note  is  also  a  combined  effort.  Capt  Abel  was 
also  responsible  for  coordinating  the  effort.  Maj  Mill  was  the  principal  author 
of  Chapter  5,  and  also  helped  with  Chapters  1,  2,  3,  4,  and  6.  Capt  Linn  wrote 
most  of  Chapters  1,  2,  3,  6,  8,  and  9.  Lt  Brown  wrote  Chapter  7. 
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Chapter  1 


INRODUCTION 

The  USAFETAC  raytrace  model  RAYTRA  is  a  computer  program  which  calculates, 
using  geometric  optics,  numerous  elements  related  to  atmospheric  refraction  ef¬ 
fects  on  a  single  beam  (ray)  of  electromagnetic  energy  while  it  travels  in  the 
Earth’s  atmosphere.  It  allows  the  user  to  completely  define  an  arbitrary  re¬ 
fracted  path  in  the  atmosphere  anywhere  from  the  Earth's  surface  to  space.  The 
model  allows  for  several  input/output  options  concerning  atmospheric  data,  path 
geometries,  and  intermediate  level  trajectory  information.  Furthermore,  output 
data  may  either  be  in  user-readable  printed  form,  machine-readable  binary  form, 
or  both. 

The  model  is  unique  concerning  its  (1)  flexibility  of  application,  and  (2) 
special  numerical  techniques  which  enable  it  to  compute  types  of  ray  paths  which 
some  models  cannot  handle.  Furthermore,  the  code  is  structured  in  a  modular, 
"top  down"  fashion  to  allow  for  ease  in  modification  and  program  maintenance.  It 
also  has  the  capability  to  produce  input  information  for  execution  of  the  Air 
Force  Geophysics  Laboratory  ( AFGL )  LOWTRAN  (low  resolution  transmittance)  pro¬ 
gram.  AFGL  TR  80-0067  [6]  describes  the  latest  versions  of  LOWTRAN  (Note:  The 
AFGL  program,  designated  at  USAFETAC  as  LOWTRN  is  very  dynamic  in  the  sense  that 
it  is  usually  in  a  state  of  constant  revision,  with  new  versions  being  issued 
periodically) . 

This  technical  note  does  not  elaborate  on  LOWTRAN  or  the  various  USAFETAC 
feeder  programs.  It  focuses  on  the  description  of  the  program  RAYTRA,  which 
performs  the  raytraces  and  contains  the  refraction  algorithms.  Also  described  in 
some  detail  are  programs  BLDATM  and  BLDCOM  which  must  be  executed  prior  to  RAYTRA 
to  set  up  input  data  and  geometric/electrical  parameters,  respectively.  Chap¬ 
ter  7  is  a  detailed  description  of  how  to  run  these  three  programs  on  the 
USAFETAC  computer  resources. 


Chapter  2 


GENERAL  DESCRIPTION  OF  RAYTRA 


2.1  Program  Descriptions 

Brief  descriptions  of  the  feeder  and  primary  programs  follow.  Figure  1  is  a 
flowchart  outlining  their  execution  sequence. 

ENAP RECON  -  processes  USAFETAC  DATSAV  upper-air  data  to  correct  many  data 
elements  [1]. 

ENAEXTR  -  Extracts  and  reformats  desired  upper-air  data  from  taped  output  of 
ENAPRECON  [2]. 

ENS2AMOD  -  Produces  USAFETAC  Point  Analysis  upper-air  data  [3]  [4]. 

PASELECT  -  Extracts  and  reformats  desired  Point  Analysis  upper-air  data  from 
output  of  ENS2AMOD  [4] . 

RKLOW  -  Optional  program  which  filters  output  of  ENAEXTR  to  selected  months, 
observation  times  (internal  documentation  only  (5]). 

BLDATM  -  Accepts  input  atmospheric  data  from  PASELECT,  ENAEXTR,  RKLOW,  disk 
files  (model  atmospheres),  or  terminal.  Processes  these  data  for  input  to  RAYTRA 
or  LOWTRN.  Calculates  and/or  models  moisture  variables. 

BLDCOM  -  Accepts  terminal  input  of  geometric  and  electrical  specifications 
and  formats  commands  for  use  in  RAYTRA  or  LOWTRN. 

RAYTRA  -  Performs  actual  raytracing  from  output  of  BLDATM  and  BLDCOM.  Con¬ 
tains  the  "raytrace  model." 

LOWTRN  (currently  LOWTRAN  5  model)  -  Performs  low  spectral  resolution  trans¬ 
mittance  and  radiance  calculations.  Accepts  input  from  BLDATM  and  BLDCOM  (6). 
FASCOD  [7]  may  be  substituted  for  LOWTRAN  for  certain  applications. 

Most  files  produced  by  these  programs  have  certain  required  names.  Detailed 
descriptions  of  the  files  and  their  formats  are  contained  in  Chapter  7. 

2.2  Assumptions  and  Limitations 

Many  assumptions/limitations  enter  into  the  raytrace  model.  Several  of  the 

more  important  ones  are 
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Figure  1.  Basic  Flowchart  of  Program  Sequence. 


a.  The  model  is  "one  dimensional"  in  the  sense  that,  for  any  single  trace, 
only  one  vertical  profile  of  atmospheric  data  is  allowed  for  input  (i.e.,  it 
assumes  spherical  or  horizontal  homogeneity  of  the  atmosphere).  Thir  is  the  most 
limiting  assumption,  although  in  many  cases  more  extensive  atmosj  c  data  is 
not  available  even  if  it  were  possible  to  use  it  as  input. 

b.  The  theory  of  geometric  optics  assumes  that,  (1)  fractior  changes  in 
the  refractive  index  (n)  within  one  wavelength  will  be  small  compar  '  1.0,  and 
(2)  fractional  changes  in  the  spacing  between  "hypothetical"  adjace.  .  rays  with¬ 
in  one  wavelength  will  be  small  compared  with  1.0  [8].  The  first  condition 
implies  that  the  refractive  index  does  not  change  appreciably  in  a  wavelength. 
The  second  condition  implies  that  ray  patterns  that  diverge  or  converge  rapidly 
or  cross  each  other  have  questionable  validity. 

c.  All  refraction  calculations  are  referenced  to  the  speed  of  light  in  a 
vacuum,  designated  as  c.  When  applying  the  model  to  compute  elements  such  as 
range  error  of  a  radar,  the  range  results  will  be  miscalculated  if  the  range  dis¬ 
play  of  the  radar  in  question  is  designed  to  electronically  compensate  for  the 
speed  of  light  in  an  atmospheric  medium  such  as  a  standard  atmosphere. 

d.  The  model  is  highly  sensitive  to  small  changes  in  refractivity  in  the 
lower  portion  of  the  atmosphere  when  computing  actual  ranges  (especially  with 
small  elevation/large  zenith  input  angles).  However,  the  range  and  angle  errors 
computed  for  a  given  refractivity  case  do  not  reflect  this  sensitivity. 

e.  The  formulas  for  computing  refractive  moduli  (N)  restrict  use  of  the 
model  to  frequencies  between  30  KHz  (wavelength  10  km)  and  1500  THz  (wavelength 
0.2  micrometers).  For  frequencies  between  115  GHz  (wavelength  0.25  centimeter) 
and  15  THz  (wavelength  20  micrometers),  the  model  results  should  be  accepted  with 
caution.  The  refractive  modulus  formulas  technically  are  not  applicable  in  that 
range,  though  the  error  may  be  quite  small.  Work  is  in  progress  at  AFGL  to  up¬ 
date  these  calculations. 


f.  The  model  assumes  that  a  ray  path  exists  between  points  of  interest. 
Tracing  terminates  if  an  apogee  point  ("ray  ducting")  or  tangent  point  (described 
below)  is  encountered,  or  the  path  intersects  the  Earth,  or  the  ray  enters  space. 
Except  in  unusual  cases,  this  signifies  that  no  path  actually  exists  between 
those  points. 

g.  Ionospheric  effects  are  not  currently  considered;  however,  provision  is 
made  for  including  an  ionospheric  module  in  future  versions.  No  refractive 
effects  are  considered  above  50  km. 

2.3  Input  Atmospheric  Data  Categories 

Input  atmospheric  data  are  grouped  into  four  main  categories 

a.  Upper-air  soundings  (i.e.,  raob  information  from  the  USAFETAC  Data  Save 
(DATSAtf)  data). 

b.  Upper-air  computations  from  the  USAFETAC/AFGWC  Point  Analysis  program 

[4]. 

c.  Special  model  atmospheres  [9]. 

d.  User-specified  upper-air  data  provided  to  the  computer  terminal  at  time 
of  execution. 

2.4  Input  Geometry  Categories 

Input  geometry  specifications  fall  into  two  main  categories 

a.  From  one  point  to  another. 

b.  From  one  point  in  the  Earth's  atmosphere  to  a  distant  celestial  object 
(one  so  far  from  the  Earth  that  rays  emitted  from  it  may  be  treated  as  parallel 
by  the  time  they  reach  the  Earth). 


There  are  two  subcategories  of  the  point-to-point  case:  upward  and  downward 
going  paths.  When  dealing  with  downward-directed  paths  from  one  point  to  anoth¬ 
er,  there  may  also  be  two  solutions:  either  a  long  path  (two  segments)  or  a 
short  path  (one  segment).  If  the  long  path  is  desired,  -ariable  flag  LEN  is  set 
equal  to  1,  for  short  paths  LEN=0.  Figure  2  depicts  the  various  geometric  cate¬ 
gories.  In  general,  two  segment  paths  are  possible  for  two  different  geometries. 
If  H2  lies  between  HI  and  the  ray  trajectory  minimum  height  H  MIN,  which  is  also 
called  the  tangent  point  (see  Figure  2c  and  2d),  level  #2  will  be  crossed  twice. 
Or  H2  may  be  specified  as  the  tangent  point  and  the  program  will  calculate  the 
trajectory  from  HI  to  H2,  and  then  from  H2  to  H3,  the  specified  ending  height. 


Figure  2.  (a)  Upward  (6^90°)  and  Downward  (0^90°)  Paths 

from  One  Point  at  HI  (HI')  to  Another  Point  at  H2  (H2’). 

(b)  Upward  and  Downward  Paths  from  a  Point  at  HI  (HI')  to  a 
Distant  Celestial  Object  (“).  Part  (c)  Shows  the  Short  Path 
Solution  for  the  Downward  Path  in  (a)  when  HMIN  §  H2  §  HI, 

Where  HMIN  is  the  Point  at  Which  Ray  is  Closest  to  Earth 
(tangent  point).  H3  is  Specified  When  Limb  Viewing  is  Desired 
(d)  and  H2  is  Set  Equal  to  the  Known  Tangent  Point.  (Extracted 
with  modifications  from  AFGL-TR-8 0-0067  [6]). 


In  addition,  the  point-to-point  cases  may  be  divided  into  those  in  which  the 
initial  ray  zenith  angle  is  either  known  or  unknown.  If  unknown,  special  treat¬ 
ment  is  required  and  is  discussed  in  Chapter  5.  Also  see  Section  4.3  for  a  more 
detailed  account  of  which  input  geometry  variable  combinations  are  acceptable. 

2 . 5  Output  Options 


Format  options  include  tape  (binary  file),  printed,  or  both.  Any  combina¬ 
tions  of  the  following  three  groups,  with  one  output  group  per  page,  is  possible 


a.  Input  data,  final  path  results,  and  informational  or  error  messages. 


b.  Intermediate  path  results  for  each  level  at  which  computations  are  made 
as  the  "ray"  is  being  traced  (at  heights  of  input  atmosphere  plus  certain  other 
critical  heights ) . 

c.  Atmospheric  data,  and  associated  refractive  moduli  (as  described  in  Chap¬ 
ter  3)  and  their  gradients. 

2.6  Input/Output  Variable  Definitions 

In  most  applications  of  the  model,  interest  will  probably  center  on  the  final 
path  results  (primarily  total  angle  error  and  range  error  when  dealing  with  radar 
systems).  The  elements  used  as  initial  input,  or  computed  in  the  final  results 
are  defined  below  and  are  depicted  in  Figure  3  for  an  upward  path  from  one  point 
to  another  and  in  Figure  4  for  a  short  downward  path  from  one  point  to  another. 

For  paths  to  distant  celestial  objects  special  definitions  of  path  results  apply 

and  are  described  later  in  this  section. 

Definitions  of  elements  in  the  final  path  results  for  tracing  from  one  point 
to  another  are  (Figures  3  and  4) 

a.  HI  -  Altitude  (MSL)  of  starting  point  of  ray,  expressed  in  kilometers. 

b.  H2  -  Altitude  (MSL)  of  ending  point  of  ray,  expressed  in  kilometers. 

c.  ANGLE  (90°-<v )  -  Zenith  angle  of  starting  point  of  ray,  expressed  in 
degrees.  Angle  a  is  the  elevation  angle  of  the  ray. 

d.  RANGE  -  Total  "apparent"  path  length  of  ray  (sum  of  computed  electrical 
retardation  along  path  length  and  actual  path  length),  expressed  in  kilometers. 
The  term  "apparent"  means  as  observed  with  a  system  such  as  a  radar. 

e.  BETA  (p)  -  Total  angle  subtended  by  path  at  the  Earth's  center,  expressed 
in  degrees.  (NOTE:  Because  the  model  traces  by  using  the  radius  of  curvature  of 
the  local  geoid,  this  angle  is  not  exact,  but  the  error  is  very  small  (see  Chap¬ 
ter  4 ) . 


f.  ARRIVAL  ANGLE  -  Zenith  angle  of  ray  when  it  reaches  the  ending  point, 
expressed  in  degrees.  Note  the  difference  in  orientation  of  this  and  the  initial 
ANGLE  (para,  c  above). 

g.  TOTAL  BENDING  (<|>)  -  Acute  angle  formed  by  intersection  of  lines  tangent 
to  ray  at  starting  and  ending  points,  expressed  in  degrees,  positive  for  normal 
curvature  (concave  downward). 
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RANGE  =  PATH  LENGTH  +  RETARD 
CORRECTION 


TOTAL  RANGE  ERROR  =  RANGE  - 
GEOMETRIC  RANGE 

GEOMETRIC  ERROR  =  FaTH  LENGTH  - 
GEOMETRIC  RANGE 


EARTH  CENTER 


Figure  3.  Upward  Path  Geometry  from  a  starting  Point  at 
Altitude  HI  (in  this  case,  HI  =  0)  To  an  Ending  Point  at 
Altitude  H2. 


h.  PATH  LENGTH  (R)  -  Length  of  curved  (refracted)  ray  path  from  starting 
point  to  ending  point  (actual  distance  energy  travels;  does  not  include  electri¬ 
cal  retardation),  expressed  in  kilometers.  NOTE:  RANGE  =  PATH  LENGTH  +  RETARD 
CORRECTION. 

i.  GEOMETRIC  RANGE  (GR)  -  Straight-line  distance  between  starting  point  and 
ending  point  of  ray,  expressed  in  kilometers. 

j.  GROUND  RANGE  (rg)  -  Great  circle  distance  along  the  Earth's  surface  from 
a  point  directly  beneath  starting  point  of  ray  to  a  point  directly  beneath  ending 
point  of  ray,  expressed  in  kilometers.  (NOTE:  The  model  treats  the  Earth  as  an 
oblate  spheroid;  therefore,  this  result  is  only  a  close  approximation;  see  dis¬ 
cussion  in  Chapter  4). 

k.  GEOMETRIC  ANGLE  (90°  -  c<G)  -  Zenith  angle  formed  by  line  representing 
GEOMETRIC  RANGE  at  the  ray  starting  point,  expressed  in  degrees.  Angle  «G  is  the 
geometric  elevation  angle  (see  Figure  8). 


V  - 


GEOMETRIC  RANGE 


Figure  4.  Downward  Short  Path  Geometry  from  a  Starting 
Point  at  Altitude  HI  To  an  Ending  Point  at  Altitude  H2 
(in  this  case,  H2  =  0). 


l.  RETARD  CORRECTION  -  "Apparent"  added  distance  ray  travels  along  PATH 
LENGTH  due  to  propagation  speed  along  the  path  being  less  than  the  speed  of  light 
in  a  vacuum. 

m.  GEOMETRIC  ERROR  (Rg)  -  PATH  LENGTH  minus  the  GEOMETRIC  RANGE. 


n.  TOTAL  RANGE  ERROR  (TRE)  -  RANGE  minus  the  GEOMETRIC  RANGE. 


o.  TOTAL  ANGLE  ERROR  (TAE)  -  Acute  angle  formed  by  the  intersection  of  the 
GEOMETRIC  RANGE  line  and  the  line  tangent  to  the  ray  at  the  starting  point,  ex¬ 
pressed  in  degrees.  Calculated  as  GEOMETRIC  ANGLE  minus  ANGLE,  it  is  usually 
referred  to  as  the  elevation  angle  error. 

p.  Earth  Radius  (RE)  -  A  constant  under  the  spherical  Earth  assumption,  but 
a  function  of  latitude  when  treating  the  Earth  as  an  oblate  spheriod. 

q.  H3  -  Altitude  (MSL)  of  ending  point  of  a  two-segment  tangent  raytrace. 
The  value  of  H3  must  be  greater  than  H2  and  flag  LEN  must  be  set  to  1  (see  Fig¬ 
ure  2 ) . 

For  paths  to  distant  celestial  objects  (those  so  far  from  the  Earth  that 
adjacent  rays  are  treated  as  being  parallel  to  each  other),  the  definitions  of 
final  path  results  and  corresponding  figures  reference  above  must  be  adjusted. 
In  fact,  the  only  final  results  of  interest  will  be  the  GEOMETRIC  ANGLE  and  the 
TOTAL  ANGLE  ERROR.  Even  though  the  other  elements,  such  as  RANGE  etc.,  are  com¬ 
puted  from  HI  to  the  "top"  of  the  atmosphere,  they  essentially  do  not  have  any 
bearing  on  this  type  of  problem.  As  such,  only  GEOMETRIC  ANGLE  and  TOTAL  ANGLE 
ERROR  will  be  redefined  below  and  depicted  in  Figure  5.  Note  that  for  type  4 


Figure  5.  Geometry  for  Path  to  Distant  Celestial  Object. 
In  this  case,  H2  must  be  the  "top"  of  the  atmosphere. 


9 


cases,  input  ANGLE  is  used  as  the  Geometric  Angle  and  the  apparent  zenith  angle 
is  calculated  using  a  Newton-Raphson  iteration  scheme. 

a.  GEOMETRIC  ANGLE  (GA^ )  -  Zenith  angle  formed  by  line  representing  GEOMET¬ 
RIC  RANGE  at  the  ray  position  HI,  expressed  in  degrees.  Note  that  this  is  the 
value  input  as  ANGLE  for  type  4  paths. 

b.  EXOGEOMETRIC  ANGLE  (GA2)  -  Zenith  angle  formed  by  intersections  of 
straight  line  (exospheric)  ray  path  and  zenith  above  HI,  expressed  in  degrees. 

c.  TOTAL  ANGLE  ERROR  (TAE)  -  Acute  angle  formed  by  the  intersection  of  the 
GEOMETRIC  RANGE  line  and  the  line  tangent  to  the  ray  at  HI.  It  is  theoretically 
equivalent  to  the  TOTAL  BENDING  term  (i|/)  defined  in  paths  from  one  point  to 
another,  expressed  in  degrees.  This  is  because  the  exospheric  ray  paths  are 
treated  as  parallel  for  an  object  so  distant  from  the  Earth. 


Chapter  3 


REFRACTION  THEORY 

3 . 1  Index  of  Refraction 

The  atmospheric  index  of  refraction  n  is  defined  as 

(1) 

where  c  is  the  electromagnetic  propagation  speed  in  a  vacuum  (2.997930  x  10  m 
sec-1)  and  v  is  the  corresponding  speed  in  the  atmosphere.  Since  the  speed  of 
propagation  in  the  atmosphere  is  slightly  smaller  than  it  is  in  a  vacuum,  n 
routinely  has  values  slightly  larger  than  unity  (near  sea  level,  approximately 
1.0003).  A  more  convenient  method  of  representing  n  is  to  express  it  as  a 
refractive  modulus  N,  where 

N  =  (n  -  1)  x  106  (2) 

At  sea  level  the  atmospheric  refractivity  ranges  from  250  to  450  N-units. 
Other  moduli  which  are  based  on  similar  types  of  transformations  also  can  be 
used.  One  such  modulus,  M,  is  related  to  N  by 

M  =  N  +  157  x  2  (3) 

where  Z  is  the  MSL  altitude  at  which  N  is  computed  and  is  expressed  in  kilome¬ 
ters.  One  advantage  of  working  with  M  instead  of  N,  especially  in  raytrace  plot¬ 
ting  programs,  is  that  the  Earth's  curved  surface  relative  to  a  refracted  beam  of 
energy  may  be  depicted  as  a  straight  horizontal  line.  For  his  model  either  N  or 
M  is  suitable  since  no  plotting  occurs.  In  fact,  the  model  uses  N. 

The  refractive  modulus  in  the  Earth's  atmosphere  decreases  approximately 
exponentially  with  height.  Since  a  "ray"  of  electromagnetic  energy,  when  passing 
through  a  medium  with  a  variable  refractive  modulus,  will  bend  toward  the  portion 
of  the  medium  with  a  higher  refractive  modulus,  a  typical  ray  trajectory  through 
the  entire  troposphere  bends  downward  toward  the  Earth.  Obviously,  the  ray  tra¬ 
jectory  through  a  smaller  segment  of  the  atmosphere  can  bend  strongly  downward, 
upward,  or  even  sideways  depending  on  the  atmospheric  refractive  modulus  gradient 
in  that  segment.  For  the  purposes  of  this  model,  the  vertical  gradient  of  N  is 
fundamental  to  the  raytracing  computations. 

3.2  Snell's  Law 

Since  this  model  assumes  that  the  refractive  modulus  in  the  troposphere 
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changes  only  with  height  above  a  spherical1  Earth  (i.e.,  horizontal  homogeneity 
is  assumed),  Snell's  Law  of  refraction  may  be  expressed  as  [10] 

n^r^  cos  =  n2r2  cos  a2  =  constant  (4) 

where  n1  is  the  index  of  refraction  at  a  reference  altitude  HI,  r^  =  HI  +  RE,  (r^ 
is  the  radius  of  the  layer  of  interest's  inner  spherical  boundry,  RE  is  the  radi¬ 
us  of  the  Earth),  is  the  ray  elevation  angle  at  HI,  n2  is  the  index  of  refrac¬ 
tion  at  some  other  altitude  H2,  r2  =  H2  +  RE,  (r2  is  the  radius  of  the  outer 
spherical  boundary),  and  a2  is  the  ray  elevation  angle  at  H2  (see  Figure  6). 


«,AT  LEVEL  2 


EARTH  CENTER 


Figure  6.  Depiction  of  Snell's  Law  Geometry  for 
Spherical  Surfaces. 


1  Actually,  the  radius  of  curvature  of  the  local  geoid,  RC,  is  used  in  the  trac¬ 
ing  equations  (see  Chapter  4).  For  purposes  of  this  discussion,  we  assume  a 
spherical  Earth. 


Since  n  (or  N)  must  be  known  to  calculate  ray  paths  from  Snell's  taw,  in¬ 
vestigators  have  derived  empirical  relationships  between  measured  atmospheric 
elements  and  N.  For  propagation  frequencies  between  30  kHz  (10  km)  and  115  GHz 
(0.25  cm),  the  well-known  Smith-Weintraub  equation  for  computing  N  from  atmos¬ 
pheric  pressure,  temperature,  and  vapor  pressure  is  generally  used.  This  equa¬ 
tion  is 

N  =  77.6  (!)  +  3.73  x  105  (%)  (5) 
i  Tz 

where  P  is  pressure  in  millibars,  e  is  vapor  pressure  in  millibars,  and  T  is  ab¬ 
solute  temperature  in  degrees  Kelvin.  Since  both  P  and  e  normally  decrease  with 
height  above  sea  level,  N  also  decreases  with  height.  In  general,  raytracing  is 
primarily  determined  by  the  gradient  of  N  rather  than  by  the  value  of  N  itself. 
For  optical  and  near-infrared  frequencies  between  about  15  THz  (20  micrometers) 
and  1500  THz  (0.2  micrometers),  an  accepted,  wavelength-dependent  equation  for  N 
is  [15] 

N  =  77.6  (|)  +  0.584  (-^~)  (6) 

1  TAZ 

where  P  and  T  are  the  same  as  in  Equation  (5)  and  A  is  the  propagation  wavelength 
in  micrometers.  There  are  no  generally  accepted  equations  for  N  between  115  GHz 
and  15  THz;  however,  the  model  assumes  Equation  (5)  is  valid  from  30  kHz  through 
3  THz  (100  micrometers)  and  Equation  (6)  is  valid  for  frequencies  fiom  3  THz 
through  1500  THz.  The  error  in  doing  this  is  not  well-known  and  great  care  must 
be  exercised  when  applying  the  model  between  115  GHz  and  15  THz.  It's  best  used 
in  the  windows  of  relative  transparency,  and  errors  of  20  N-units  are  common  even 
then. 

Equations  (5)  and  (6)  are  used,  given  input  atmospheric  data,  to  compute  N 
for  each  MSL  height  associated  with  the  atmospheric  data.  For  the  vertical 
legion  between  the  uppermost  level  of  atmospheric  data  and  50  km,  N  is  assumed  to 
decay  exponentially  from  the  last  computed  value  to  a  value  of  zero  at  50-km  al¬ 
titude.  It  is  assumed  that  no  refraction  occurs  above  50  km  (i.e.,  ionospheric 
effects  are  not  currently  considered).  Above  50  km,  data  levels  are  set  at 
increments  of  DH,  which  has  a  default  value  of  50  km.  At  the  DH  default  value  of 
50  km  and  using  a  32-level  atmosphere  up  to  the  50-km  level,  the  prog tain  is  only 
set  up  to  handle  H2  or  H3  values  of  4000  km  or  less,  because  of  data  array  limi¬ 
tations.  Of  course,  above  50  km  the  raytracing  consists  only  of  a  stiaight  line. 
If  larger  H2  and  H3  values  are  input,  the  value  of  DH  must  be  increased  accord- 


Chapter  4 


BASIC  GEOMETRIC  CONSIDERATIONS 


This  chapter  describes  the  mathematical  basis  for  the  ray  tracing  calcula¬ 
tions.  The  equations  are  primarily  based  on  the  geometry  of  the  tropospheric 
propagation  problem  and  small  number  approximations.  Some  of  the  numerical  tech¬ 
niques  used  in  the  model  are  described  in  Chapter  5.  The  reader  should  acquire 
from  Chapter  4  an  understanding  of  the  general  framework  about  which  the  details 
of  the  model  are  constructed.  Some  form  of  the  equations  are  also  found  in 
references  [8],  [10],  and  [12]. 


Basic  Tracing 


lations 


Consider  a  spherical  Earth  with  a  reference  level  just  above  the  surface  of 
radius  r^  =  RE  +  HI.  Consider  a  spherically  concentric  tropospheric  layer  above 
the  Earth  with  a  top  at  H2  and  a  radius 


where 


r2  =  RE  +  H2  =  rL  +  AH 
AH  =  H2  -  HI 


and  a  ray  starting  at  HI  to  an  ending  point  at  H2  along  a  typical  refracted  path. 
Further,  consider  that  the  vertical  gradient  of  the  refractive  modulus  is  con¬ 
stant  in  the  layer,  and  the  ray  has  a  starting  point  elevation  angle  of  c^,  an 
ending  point  elevation  angle  of  a2  and  a  straight-line  distance  (GEOMETRIC  RANGE) 
between  HI  and  H2  of  GR  (see  Figures  7  and  6).  Snell's  Law  may  be  rewritten  from 
Equation  (4)  as 


njr^  cos  ai  =  n2  ^ri  +  AH)  cos  a2  (7) 

or 

(1  +  N1  x  10-6)  r^  cos  =  (1  +  N2  x  10-&)  (r^  +  AH)  cos  c»2  (8) 

If  a  and  b  are  small  quantities  compared  to  1.0  the  approximations 

(1  ±  a)m  ~  1  ±  ma  (9) 


and 

(1  ±  a)m  (1  ±  b)n  ~  lfmalnb 


(10) 


may  be  used.  Remember  that  N  x  10"6  is  on  the  order  of  0.0003  and  even  for  an 
extremely  large  AH  value  of  50  km,  AH  *  ^  is  approximately  0.008.  Using  Equa¬ 
tion  (10)  we  have 
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Figure  7.  Depiction  of  Basic  Geometry  for  a  Refracted  Path 
Between  Two  Points  Through  One  Spherical  Layer.  Also  included 
is  the  derivation  of  Equation  (15). 


(1  +  N2  x  10"6)  r  (1  -  Nx  x  10-6)  :  1  -  (Hj  x  10"6)  +  (N2  x  l(Tb)  (11) 
Simple  algebra  can  be  used  to  write 

(rx)  =  (r1  +  AH)  =  [1+(AH  *  rx ) ]  _1  (12) 

Using  the  relationship  in  Equations  (10)  and  (11),  Equation  (8)  may  be  rewritten 


cos  a2  ~  (cos  a^)  r  [(1  +  AH  t  r^)  (1  -  (N^  -  N2 )  x  10 


or  simplied  further  to 


cos  u2  ~  (1  +  (N^  -  N2)  x  10~  -  (AH  r  r^ ) |  (cos  1 1^)  (14) 

Note  that  a  is  dependent  only  on  height  AH  (which  is  the  distance  between  the 
levels  of  interest)  and  (N^-N2),  which  is  the  difference  in  the  refractive  indi¬ 
ces  between  the  HI  and  H2  levels  of  interest.  Thus,  n  is  independent  ol  the 
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shape  of  the  refractive  profile  between  HI  and  H2.  This  fact  allows  us  to  deter¬ 
mine  a  at  each  level  of  interest  independently  of  other  levels  [8]. 

It  is  clear  from  Figure  7  that  the  elevation  angle  a  is  related  to  the  total 
bending  tji  and  the  Earth-centered  angle  p.  Simple  geometry  gives  the  relationship 

p  =  «|i  +  a2~ai  (15) 

An  equally  important  relationship,  derived  by  Weisbrod  and  Anderson  |10],  is 
the  dependence  of  the  amount  of  ray  bending  on  the  vertical  refractive  modulus 
gradient 


=  [2  (N^  -  N2)  x  10-6)  t  (tan  +  tan  a2 ) 


(16) 


where  i|<  is  expressed  in  radians.  Total  bending  is  the  sum  of  the  layer  contribu¬ 
tions  . 

Refractive  effects  of  the  atmosphere  cause  propagation  delays  and  range 
errors.  First  consider  the  geometric  error  (Rg)  defined  as  the  difference  be¬ 
tween  the  actual  path  length  (R)  and  the  geometric  range  (GR)  (see  Figures  3  and 
4).  It  can  be  shown  geometrically  using  the  law  of  cosines  that 

GR  =  [ (RE+H1 )2  +  ( RE+H2 ) 2  -  2(RE+H1)  ( RE+H2 )  cos  p]*5  (17) 

Distance  R  is  found  similarly  by  using  the  form  of  Equation  (17)  and  summing 
incremental  calculations  using  incremental  changes  in  height  that  are  smaller 
than  H1-H2  =  AH  and  corresponding  smaller  p  values.  If  the  N  gradient  is  con¬ 
stant  between  HI  and  H2,  then  GR  =  R.  In  summary  we  may  write 

Rg  =  R  -  GR  (18) 

Another  refractive  effect  experienced  by  a  beam  of  electromagnetic  energy  in 
the  atmosphere  is  electrical  retardation.  This  represents  the  slowing  of  propa¬ 
gation  along  the  refracted  path  and  is  independent  of  bending.  If,  for  example, 
a  beam  of  energy  travels  through  a  homogeneous  atmosphere  in  which  the  gradient 
of  N  is  zero,  the  beam  will,  because  of  molecular  interaction,  travel  slower  than 
it  would  in  a  vacuum.  From  Equation  (1),  we  have,  given  v  =  ds/dt, 

(1  +  N  x  10-6)  ds  =  c  dt  (19) 

From  Equation  (19)  one  can  see  that,  when  c  is  employed' to  compute  the  distance 
traveled  by  propagating  energy  in  a  given  amount  of  time,  an  error  of  N  times  106 
times  true  range  will  result.  For  a  range  path  increment,  Rl,  defined  by  levels 
1  and  2,  with  a  constant  N  gradient,  the  retardation  error  is 


16 


and  the  summation  over  the  entire  path  of  the  incremental  error  valuer  will  be 
called  Re.  Consequently,  the  total  range  error,  TRE,  will  be  the  sum  of  Re  and 
the  added  geometric  distance,  Rg,  along  the  total  refracted  path  length 


TRE  =  Re  +  Rg  (21) 

Typically  the  value  of  Re  is  much  larger  than  the  value  of  Rg.  More  detailed 
discussions  of  these  basic  concepts  are  given  in  references  [8],  ( 10 1 ,  and  1 12  J . 

Two  additional  geometric  relationships  are  depicted  in  Figure  8  |9|  and  may  be 
written  as  follows 


and 


cos 


-1  [ (RE  +  HI)  -  (RE  +  H2)  COS  8  j 


GR 


(22) 


AE  =  cr^  -  aG  ( 23  ) 

where  aQ  is  the  elevation  angle  of  a  nonrefracted  path  GR,  and  AE  is  the  eleva¬ 
tion  angle  error  (all  angles  expressed  in  radians). 


Figure  8.  Depiction  of  the  GEOMETRIC  RANGE,  GR,  and  Its 
Relationship  to  the  Elevation  Angle  Error,  AE. 


Figures  6,  7,  and  8,  depict  the  simple  geometry  for  bending  through  one  atmo¬ 
spheric  layer  concentric  to  a  spherical  Earth.  In  practice,  the  ray  must  tra¬ 
verse  many  successive  layers  in  the  atmosphere  defined  by  reported  levels  of 
atmospheric  data.  In  addition,  the  model  treats  the  Earth  as  an  oblate  spheroid. 
The  consequences  of  this  are  outlined  in  the  following  paragraphs. 

4.2  Earth  Radius  and  Radius  of  Curvature 


The  form  of  Snell's  Law  used  by  the  model  (Equation  (7))  assumes  a  spherical 
Earth,  but  the  critical  point  is  not  the  Earth  radius  per  se,  but  the  curvature 
of  the  surface .  Therefore,  the  radius  required  by  Equation  (7)  is  the  radius  of 
curvature  of  the  Earth's  surface  in  the  vicinity  of  the  ray  path.  To  a  second 
approximation,  the  Earth  is  an  oblate  spheroid  and  RAYTRA  assumes  this  in  various 
calculations  (Figure  9). 


MOUTH  POLE 
I 


I 


I 


EQUATOR 


Figure  9.  Depiction  of  the  Earth  as  an  Oblate  Spheriod. 

The  symbol  0  is  the  latitude  at  point  P,  RE  is  the  computed 
Earth  radius  at  point  P,  and  RC  is  the  radius  of  curvature 
of  the  local  geoid  about  point  P. 


RAYTRA  calculates  the  radius,  RE,  and  the  radius  of  curvature,  RC,  for  the 
latitude  of  the  atmospheric  data  or,  if  the  latitude  is  unknown,  a  default  lati¬ 
tude  of  45  degrees.  The  radius  is  computed  by  [13] 

RE  =  A(1  -  f  sin2  <t>)  (24) 

where  A  is  the  semimajor  (equatorial)  radius  (6378.139  km)  and  f  is  a  constant., 

approximately  equal  to  (A  -  B)/A,  with  a  currently  accepted  value  of  3.35282  x 
-3  .  . 

10  ;  B  is  the  semiminor  (polar)  radius  (6356.750  km)  and  <(>  is  the  latitude. 

The  radius  of  curvature  of  an  oblate  spheroid  is  given  by  [ 12 | 
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? 

RC  - _ A _ _ _ _  (2b) 

B  {1  t  C  COS2  (f)"5  (1  +  C  COS2  <t>  cos2  0) 

2  2 

where  A,  B,  and  0  are  defined  above;  C  =  A  /B  -  1;  and  fl  is  the  ray  azimuth, 
clockwise  from  north.  The  model  currently  assumes  an  azimuth  of  zero  (north) 
unless  the  user  specifies  the  ray  azimuth. 

The  path  subtense,  BETA,  input  by  the  user,  is  measured  from  the  center  of 
the  Earth  and  would  normally  be  calculated  from  a  knowledge  of  the  latitudes  and 
longitudes  of  HI  and  H2 .  Since  the  tracing  algorithms  use  the  radius  of  curva¬ 
ture,  RC,  BETA  must  be  converted  on  input  and  reconverted  on  output.  The  equa¬ 
tions  used  for  these  conversions  are  discussed  below  and  depicted  in  Kiquie  10. 

If  a  value  of  BETA  is  supplied  by  the  user,  it  is  assumed  to  subtend  the 
ground  range,  rg,  from  the  center  of  the  Earth  (p^  ).  The  ground  range  rg  is 
held  constant,  and  the  subtense  from  the  center  of  curvature  (p  )  is  calculated 
by 

fc  '  Pin  <M>  1261 

Once  tracing  is  complete,  the  subtense  must  be  reconverted.  It  is  assumed 
that  the  user  is  interested  in  the  value  from  HI  to  H2,  shown  as  p  (  in  Fig¬ 
ure  10.  From  the  Law  of  Cosines,  this  is 


cos  <POUt> 


[RE  +  HI)2  (RE  +  H2 ) 2 
2 {RE  +  HI)  (RE  +  H2 ) 


This  choice  for  the  output  value  of  BETA  is  somewhat  arbitrary,  but  seems  the 
best  choice;  the  difference  between  it  and  p.  is  very  small  in  any  case  (note 
that  Figure  10  is  greatly  exaggerated),  and  the  use  of  Pc  in  the  tracing  equa¬ 
tions  increases  the  overall  accuracy  of  the  results. 

Equations  ( 1 )  —  ( 6 )  in  Chapter  3  and  (7)-(27)  in  this  chapter  form  the  basis 
for  solving  the  many  different  types  of  raytracing  problems  of  which  the  model  is 
capable.  Special  iterative  techniques  are  used  in  conjunction  with  the  above- 
mentioned  equations  to  solve  those  problems  where  the  initial  zenith  angle  is 
unknown  (see  Chapter  5).  The  particular  combinations  and  sequence  ot  equations 
used  depend  on  which  input  variables  are  known  and  provided  by  the  user . 


4.3  Input  Geometry  Variables  and  Path  Types 


I 


There  are  many  possible  combinations  of  input  geometry  variables  tot  which 
the  model  will  calculate  the  refracted  path.  In  general,  those  variables  not 
supplied  (or  not  used  as  inputs  by  the  model)  are  calculated,  and  the  output  will 
consist  of  a  consistent  set  of  all  geometry  variables  as  defined  in  Section  2.6. 
These  various  combinations  are  discussed  briefly  in  the  following  paragi uphs  and 
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Figure  10.  The  Relationship  Among  the  Input,  Output,  and 
Internal  Values  of  the  Variables  BETA:  The  "Earth-centered" 

Substense  of  the  path.  The  difference  stems  from  treating 
Earth  as  an  oblate  spheroid. 

the  methods  of  calculation  outlined  more  fully  in  Chapter  5.  For  convenience  and 
ease  of  understanding,  they  will  be  grouped  by  general  methods  of  calculation. 

There  are  three  basic  path  types  used  in  RAYTRA,  which  are  selected  by  set¬ 
ting  the  variable  I  TYPE  to  2  or  3  or  4.  This  notation  is  consistent  with  that 
used  in  the  AFGL  program  LOWTRAN .  Path  type  1  is  not  used  by  RAYTRA.  Type  2  is 
a  path  between  two  points  and  encompasses  most  of  the  paths  usually  encountered 
in  refraction  problems.  Types  3  and  4  involve  paths  to  or  from  a  distant  (at 
"infinity"),  usually  celestial  object,  where  rays  arriving  at  the  "top"  of  the 
atmosphere  (nominally  50  km)  are  assumed  to  be  parallel.  Type  3  paths  require 
that  HI  and  the  apparent  (observed)  ANGLE  at  HI  be  input.  Type  4  involves  the 
true  (unrefracted)  or  ephemeris  angle  at  HI.  Otherwise,  the  geometry  is  like 
I  TYPE  3.  In  I TYPE  3  cases,  geomtric  angle  and  total  angle  error  are  calculated. 
In  I TYPE  4  cases,  apparent  zenith  angle  and  total  angle  error  are  calculated. 

Paths  of  I  TYPE  2  are  further  divided  into  those  for  which  the  ray  zenith 
angle  at  HI  is  known,  and  those  for  which  it  is  not.  Acceptable  combinations  of 


20 


input  geometry  variables  for  Type  2  paths  are  outlined  below.  Note  that  Ill  must 
always  be  supplied. 

a.  HI,  ANGLE  and  at  least  one  of  H2,  RANGE,  or  BETA,  are  supplied.  If  more 
than  one  of  the  latter  are  supplied,  all  but  one  are  ignored,  with  the  priority 
of  selection  being  the  order  given  above.  For  example,  if  HI,  ANGLE,  H2,  and 
BETA  are  given,  BETA  is  ignored.  For  all  three  possible  combinations  described 
in  this  paragraph  tracing  begins  at  HI,  at  zenith  angle  ANGLE,  and  proceeds  until 
the  condition  specified  by  the  third  variable  is  satisfied,  or  the  ray  either 
leaves  the  defined  atmosphere,  intersects  the  surface,  encounters  a  duct,  or 
reaches  a  tangent  point  (see  Figure  2).  I  f  H2  is  supplied,  tracing  terminates 
when  the  ray  reaches  an  altitute  of  H2.  If  ANGLE  is  greater  than  90  degrees 
(downward  path)  and  H2  is  less  than  HI,  tracing  terminates  on  the  first  or  second 
crossing  of  the  level  H2,  depending  on  whether  the  input  variable  LEN  (used  to 
define  longest  of  two  paths)  is  set  to  0  or  1,  respectively  (see  Figure  2).  If 
RANGE  or  BETA  are  supplied  instead  of  H2,  tracing  terminates  when  the  appropriate 
value  of  apparent  ("radar")  range  or  Earth-centered  subtense  is  reached.  This 
normally  requires  iteration  on  the  last  path  segment  (between  two  atmospheric 
data  levels)  as  outlined  in  Chapter  5. 

b.  HI  and  H2  are  supplied,  but  not  ANGLE.  The  path  is  assumed  tangent  at  H2 
(a2=0).  If  H2  is  greater  than  or  equal  to  HI,  the  program  cannot  compute  a  path 
and  an  error  message  is  written.  If  H2<H1,  the  roles  of  HI  and  H2  are  temporar¬ 
ily  exchanged  and  tracing  proceeds  from  H2  to  HI  with  an  initial  zenith  angle  at 
90  degrees.  Results  are  stored  and  written  in  the  normal  order  from  HI  to  H2. 
If  H3  has  also  been  supplied  (see  Figure  2)  and  is  greater  than  H2  and  flag 
LEN=1,  tracing  continues  beyond  H2  and  terminates  at  H3.  If  LEN  is  set  to  1  and 
H3  is  not  provided,  a  default  value  of  50  km  is  used  for  H3  and  the  long  path  is 
still  computed.  This  allows  an  arbitrary  tangent  path  to  be  defined  by  the  end 
points  HI  and  H3  and  the  tangent  height  (lowest  point  of  the  path)  H2 .  This  is 
the  only  case  in  which  the  variable  H3  has  any  meaning  to  the  tracing  al gorithms. 

c.  HI,  H2 ,  and  RANGE  are  supplied,  but  not  ANGLE.  An  initial  estimate  of 
ANGLE  is  made,  assuming  a  standard  refraction  model  (the  4/3  Earth  approximation 
[14])  and  tracing  proceeds  to  H2  as  in  paragraph  a,  above.  The  resulting  value 
of  the  apparent  range  is  compared  to  the  input,  RANGE,  and  tracing  proceeds 
iterativly  until  they  match  with  the  limit  specified  by  the  internal  program 
variable  RNGACC  (currently  0.0001  km).  The  iterative  scheme  is  discussed  in 
detail  in  Chapter  5.  To  avoid  possible  confusion  over  the  value  of  flag  LEN,  the 
input  heights  may  be  temporarily  exchanged  to  assure  that  tracing  proceeds  from 
the  lower  altitude  to  the  higher. 

d.  HI,  H2,  and  BETA  are  supplied,  but  not  ANGLE  or  RANGE.  The  procedure  is 
identical  to  that  in  paragraph  c  with  BETA  assuming  the  role  of  RANGE.  The  pre¬ 
cision  criteria  is  in  the  variable  BETACC  ( =RNGACC/RC20 . 000001  radians). 


Input  altitudes  are  instrumental  in  defining  the  "top"  of  the  atmosphere. 
Subroutine  REFMOD  will  construct  a  refractive  modulus  profile  from  the  surface 
(assumed  to  be  the  lowest  level  of  the  atmospheric  data)  to  a  height  equal  to  the 
maximum  of  HI,  H2,  H3,  or  50  km.  All  refractive  moduli  above  50  km  are  currently 
set  to  zero.  To  avoid  array  overflow,  if  the  maximum  heights  exceeds  4000  km, 
the  value  of  DH  must  be  specified  as  a  number  larger  than  the  default  of  50  km  as 
explained  in  Chapter  3 . 


Chapter  5 


SPECIAL  NUMERICAL  TECHNIQUES 


5 . 1  Cases  where  ANGLE  is  Known 


The  tracing  equations  introduced  in  Chapter  4  are  applied  in  a  straight-for¬ 
ward  manner  if  the  initial  zenith  angle  is  known.  The  elevation  angle  at  the 

next  layer  boundary  (ab)  is  calculated  by  Equation  (13).  If  cos(«b)  is  greater 

than  one,  a  ray  maximum  or  minimum  exists  somewhere  at  an  unknown  level  HM  within 
the  layer  and  is  found  by  setting  cos(ab)  =  1  and  solving  for  HM  in  place  of  H2 
in  Equation  (13)  as  illustrated  below. 

The  general  form  of  Equation  (13)  is 

cosc»a 

cos“b  =  TrT~Kz7i~YTTTTn7Wz~Kz) 

a 

where  a  =  elevation  angle 

subscript  a  -  refers  to  the  known  (previous  level)  value, 

subscript  b  -  refers  to  the  unknown  (next  level)  value. 

AZ  =  in  general  the  layer  thickness,  Hb-Ha,  in  our  case  HM-Ha 
r  =  the  distance  from  the  center  of  curvature  to  H,  (i.e.  R  +  HJ 

=  the  vertical  refractive  index  gradient;  (nb-na)/(Hb~Ha ^ 

At  a  maximum  or  minimum  point,  cosab  =  1,  giving 

(1  +¥)  (1  +  |f)  =  cosaa 

a 

rearranging 


/ra>  <AZ>2  +  +  If)  AZ  +  -  cos"a>  =  0 

a 

This  equation  may  be  written  in  the  form  of  the  quadratic  formula, 


AZ 


_  -B  ±  l  B2  -  AC 
A 


where 


_  1_  an 

r.  dz 
a 


B  =  (1/2)  (5-  +  §*> 
d 

C  =  1  -  cosa 

Q 

Of  the  two  possible  solutions,  we  always  want  the  smaller  in  absolute  magnitude. 
Rather  than  finding  both  and  testing,  we  force  calculation  of  the  smaller  by 


A  Z  =  ± 


^  =  ±  FACTOR  =  HM  -  Ha 


where  the  plus  sign  is  for  upward  paths  (i.e.,  a  maximum  point).  Rewriting,  we 
get 

HM  =  Ha  ±  FACTOR  (28) 

This  procedure  assumes  a  constant  refractive  gradient  between  H  and  H.  .  This 
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assumption  is  also  used  to  estimate  a  value  for  N  at  HM. 

Once  of ^  is  determined  (either  zero  as  discussed  above  or  some  larger  number), 
the  bending  (41)  is  calculated  by  Equation  (16),  the  subtense,  0,  by  Equation 
(15),  and  the  range  increment  within  the  layer  is  found  using  Equation  (17).  The 
retard  correction  is  found  by  Equation  (20)  the  incremental  values  for  the  layer 

are  added  to  running  sums,  and  the  procedure  is  repeated  for  the  next  layer  un¬ 

less  one  of  the  stopping  criteria  is  met.  If  the  trace  is  to  be  terminated  by  an 
altitude  (usually  H2),  the  summation  stops  when  that  point  is  reached.  If,  how¬ 
ever,  RANGE  or  BETA  are  used  as  the  stopping  criteria,  an  iteration  must  be 

performed  between  the  last  two  given  levels  bounding  the  last  layer,  since  the 
next  height,  H^',  is  unknown. 

If  the  apparent  range  is  used  to  terminate  the  trace,  the  last  between-level 
increment  is  approximated  by  a  right  triangle  of  sides  DR,  DH,  and  a  horizontal 
range  increment  (see  Figure  11a).  The  triangle  is  solved  for  DH,  the  new  layer 
is  retraced  and  the  resultant  range  increment  compared  to  the  required  increment, 
and  the  procedure  is  repeated  (Figure  lib)  until  the  error  is  less  than  a  prede¬ 
termined  value  as  given  by  the  variable  RNGACC  (currently  0.001  km).  The  proce¬ 
dure  is  analagous  for  downward  paths.  If  the  estimated  H2  is  initially  closer  to 
than  Ha,  the  iteration  starts  from  H^,  thereby  reducing  the  number  of  itera¬ 
tions  required.  If  the  ray  is  nearly  horizontal,  numerical  instability  may 
develop  and  the  calculated  step,  DH,  is  multiplied  by  a  relaxation  factor,  cur¬ 
rently  0.9. 


If  the  subtense,  0  is  used  to  terminate  the  trace,  iteration  on  the  last 
layer  proceeds  by  a  Newton-Raphson  technique.  We  may  describe  the  relationship 
between  0  and  the  height,  expressed  as  the  total  distance  from  the  center  of 
curvature,  r,  in  the  general  functional  form 


Figure  11.  (a)  The  Right  Triangle  Approximation,  Used 

Iteratively  to  Find  the  Value  of  H2  Which  Gives  the  Required 
Input  RANGE,  (b)  Second  Iteration  of  the  Right  Triangle 
Approximation  Shown  for  an  Upward  Path. 


P  =  f(r) 

Expanding  in  a  Taylor  series  about  some  rQ  (initially  rg)  and  truncating  higher 
order  terms  gives 


P  =  P0  +  If  (r  -  r0) 


Solving  for  the  height  increment,  Ar  from  level  a  gives 


Ar  =  (p1  -  pQ)  t  (||)  (29) 

where  p'  is  the  required  value  of  subtense  as  found  by  Equation  (26)  in  Chap¬ 
ter  4.  This  equation  is  applied  iteratively,  using  the  new  rQ  and  updating  pQ  by 
retracing,  until  convergence  to  within  BETACC.  The  key  to  success  is  the  calcu¬ 
lation  of  the  partial  derivative.  In  many  Newton-Raphson  procedures,  it  is  esti¬ 
mated  numerically  by  finite  differencing.  This  practice  often  results  in  slow 
convergence  and  numerical  instability.  For  this  model,  it  was  determined  analyt¬ 
ically  and  is  evaluated  by  the  program.  It  is  derived  as  follows.  We  differen¬ 
tiate  Equation  (15)  by  r,  giving 


ajs  _  a“b  ,  at*  a"a 
Sr  ~  Sr"  +  Sr  ”  SF" 


(30) 
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The  last  term  is  zero  since  a,  does  not  change.  The  first  term  on  the  right  side 

a 

is  found  by  differentiating  Equation  (7).  The  left  side  of  the  result  is  zero 
since  the  value  at  level  a  does  not  change.  This  gives  (dropping  subscript  b) 

0  =  -nr  sin  a  +  r  cos  a  +  n  cos  a 

Since  3r/3r  =  1  and  3n/3r  is  constant  within  a  layer,  we  obtain,  substituting  yn 
for  3n/3r  and  rearranging 

3a  _  cos  c(rYn  +  n) 

3r  nr  sin  a 

Dividing,  expanding,  and  discarding  small  terms 


~  (cot  ot)(i  +  yn) 


(31) 


The  second  term  on  the  right  side  of  Equation  (30)  is  found  by  differentiating 
Equation  (16)  with  An  substituted  for  (NQ  -  N)  10-6 


3¥ 


3r  tan  +  tan  a b  3r 


3(An)  + 


2(An) 


(tan  <*a  +  tan  a^) 


3 (tan  +  tan  a^) 


3r 


Again,  since  the  variables  at  level  a  are  constant,  this  simplifies  to  (dropping 
the  subscript  b) 


3Y 


-2  Yr 


3  tan  or 


3r  tan  a  +  tan  a  tan  a  +  tan  or  3r 

a  d 


-1 


tan  a,  +  tan  a  ^  ^n  + 


3a . 


sec2  a  ^) 


(32) 


where  3a/3r  is  calculated  by  Equation  (31).  Equations  (30),  (31),  and  (32)  allow 
the  iterative  solution  by  Equation  (29).  The  procedure  is  quite  stable  and  con¬ 
verges  within  three  to  five  iterations.  The  only  problem  occurs  when  tan  a  =  0 
(see  Equation  31).  In  that  case,  Ar  is  reduced  by  1  percent  and  the  procedure  is 
repeated. 


5.2  Cases  where  ANGLE  is  Unknown 


When  the  initial  zenith  angle  is  unknown,  cases  (c)  and  (d)  in  Chapter  4,  the 
general  procedure  is  to  guess  at  a  value  of  ANGLE,  perform  the  trace,  compare  the 
final  values,  compare  the  value  of  RANGE  or  BETA  with  the  desired  values,  adjust 
the  value  of  ANGLE,  and  repeat  until  within  RNGACC  or  BETACC  of  the  desired 


rr* 


values  of  RANGE  or  BETA,  respectively.  The  method  of  updating  involves  a  Newton- 
Raphson  scheme  similar  to  that  outlined  in  Section  5.1. 

5.2.1  Iteration  to  BETA.  The  truncated  Taylor  expansion  gives 


P 


where  is  the  estimated  elevation  angle  at  HI,  0e  is  the  result  of  a  trace 
using  aQ,  and  aQ'  is  the  new  guess  for  aQ.  Rearranging 


(P  -  Pe> 

“o'  “  “o  +  (3p/3aQ)  (33) 

The  calculation  of  3p/3a0  proceeds  in  manner  similar  to  that  in  Section  5.1, 
except  that  p  is  treated  as  the  sum  of  incremental  values  over  all  atmospheric 
layers  in  the  path 


n 

I 

i=l 


where  n  is  the  total  number  of  layers,  i  runs  consecutively  from  HI  to  H2  and  i  = 
1  at  HI  (aQ  =  ) .  We  will  develop  the  solution  for  one  layer,  from  i  to  i  +  1, 
which  is  evaluated  by  the  program  as  tracing  proceeds,  accumulating  a  sum  in  the 
process.  We  will  see  that  the  minimum  points  require  special  treatment  and  we 
here  define  the  index  of  the  minimum  point  as  i  =  m. 


Again,  given  Equation  (15),  we  have 


9“i+l 

3a„ 


(34) 


From  Equation  (7) 


3“0  •  9“i 

-n0r0  sin  aQ  g^  =  ^r.  sin  a.  gj- 


+  ni  cos  “i  5iJ  +  ri  cos  “i  5^ 

where  g-jj—  =  1,  of  course. 

For  i  *  m,  the  last  two  terms  are  zero  and  for  i  =  m  (the  layer  past  the  minimum 
point),  the  first  term  on  the  right  is  zero  since  a  s  0.  Therefore,  for  i  t  m, 
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Again,  from  Equation  (7) 


3o,i  _  nQrQ  sin  aQ 
3aQ  -  n^r^  sin  a^ 


n0r0  =  niri  cos  a. 


3a  ^  tan 

3aQ  ~  tan  a^ 


for  i  t  m 


Similarly 


3“i+l 


tan  a. 


3a0  tan  a.+1 


for  i  +  1  *  m 


for  i  =  m,  we  are  left  with 

noro  sin  “0  =  -ni  *  ri  ^  1  =  m  <37) 

which  will  be  used  below.  Now,  for  the  second  term  in  Equation  (34);  from  Equa¬ 
tion  (16),  we  get 


34^  3  (nt  -  ni+1)/3aQ 

3aQ  "  1/2  (tan  a^  +  tan  «^+1) 


3  [1/2  (tan  a^  +  tan  a^+1) 


1/2  (tan  a^  +  tan  a^+1) 


For  i  *  m  and  i  +  1  t  m,  the  first  term  is  zero.  Differentiating  the  tangent 
functions  and  using  Equations  (35)  and  (36),  we  get 

34>i  -<l<4  tan  an  .  , 

3«0  tan  ai  +  tan  a^+1  'cos  a^  sin  a^  cos  a^+1  sin  a^  +  1'  ’ 


for  i  *  m,  i+1  *  m. 


For  i  *  l  =  m,  (the  layer  before  the  tangent  point),  we  get 


-2  9ni  +  l  +i  tan  a0 


3a 0  tan  a^  3a  ^ 


i+1  =  m 
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By  the  chain  rule 


3n-  3n.  3r-  3r. 

_ i.  =  1 _ _  =  v  — -  (40) 

a«i  a«0  *  aa0 

where  y  is  the  gradient  from  i  to  i  +  1.  Substituting  into  Equation  (37)  and 
rearranging 


9ri 

-noro  sin  “o  =  a^  (ni  +  ri  *> 


(41) 


Back  substituting  into  Equation  (40)  gives 


9ni  _  -no  r0  sin  “0 


a« 


o  +  ri 


Substituting  into  Equation  (39)  results  in 


a*i>  d 

3a, 


2  n0r0  sin  “o 


tan  a.  (n.+1/Y  +  r.+1) 


4.i  tan  aQ 
sin2  a- 


i+1  =  m 


(42) 


For  i  =  m,  the  derivation  is  similar;  however,  since  the  atmosphere  is  spheri¬ 
cally  symmetric,  a^(i+l  =  m)  =  -  =  m)  and  the  derivative  is  symmetric 
around  the  minimum  point,  i.e., 


3a7 


i+1  =  m  = 


da. 


Consequently 


8>i 


tan  a, 


tan  a, 


30^  tan  tan  ot^+^ 


«l»i  tan  «0 


tan  a ^  +  tan  u^+1  'cos  sin 


cos  «i+1  sin  .»i  +  i 


(43) 


i  *  m,  i+1  *  m 


da 


tan  a. 


2  nQr0  sin  aQ 


0  tan  ai 


tan 


("i+l/Y  +  ri+i > 


il«i  tan  «0 


sin' 


Of  . 

1 


i  +  1  =  m 


(44) 


a*i 


tan  a 0  2  nQr0  sin  aQ  tan  ifQ 

9o"g  “  tan  or i  +  1  +  tan  «i  +  1  (n^y  +  ri )  "  sin2  <»iu 


(45) 
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5.2.2  Iteration  to  RANGE.  The  Newton-Raphson  equation  is 


«0'  =  «0  +  (R  “  Re)/(3R/3a0) 


(46) 


where  R  is  the  desired  apparent  range  and  Rg  is  the  estimate  of  R  calculated  from 
simple  geometry.  We  begin  by  rewriting  Equation  (17)  for  a  single  layer  as 

Ri  =  <ri2  +  W  -  2  riri+i  cos  Pi>H  <47> 


Differentiating  by  an  gives 


3R 


‘  ‘  2r,  “-1 


SB- 


3aQ  2Ri  i  3aQ  “i+1  3aQ 


+  2ri  ri+i  sin  h  3^7 


2ri+l  cos  Pi  3^  -  2ri  cos  Pi  ^7 


(48) 


For  i  #  m,  i+1  *  m; 


*Ji 

3a,. 


0  and 


aRi  i  3Pi 

3^7  =  S7  (ri  ri+i  Sln  Pi  55T>  1  *  m'  1+1  *  m 


(49) 


where  dp^/daQ  is  given  by  Equation  (43). 


For  i  =  m  and  i  +  1  =  m,  there  is  an  additional  term,  G(aQ). 
we  have  ( ignoring  the  term  in  Equation  49 ) 


From  Equation  (43), 


G(V  =  k~  <ri+i  -  ri cos  Pi  1  +  1  =  m 


Using  Equation  (41)  and  rearranging  gives 


G(«o>  = 


i  noro  sin  a0 

-  6-  <ri  +  l  -  ri  cos  Pi>  nr  +  r."  1  +  1  =  m  (50) 

Ri  11  1  1  'ni+l  ri+lY  v) 


Since  the  atmosphere  is  symetric,  it  can  be  easily  shown  the 


Therefore; 


C<°0>i=m  =  c<a0>i+l=ro 


3R.  .  3^ 

3^-  =  R"  <ri  ri+l  Sin  Pi  8^>  +  G(“0> 


(51) 


B 


where  G(uQ)  ~ 


i  (Eq.  50);  1  =  m,  1  +  1  =  m 
{  0  ;  i  *  m,  i  +  1  *  m 


5 . 3  Path  Types  3  and  4 

For  paths  to  distant  celestial  objects,  the  rays  beyond  modeled  atmospheric 
influences  are  assumed  to  be  parallel  (see  Figure  12). 


H2  / 

"  /" 


PARALLEL 

RAYS 


--  —  EtlTN  CENT!*  J  I 

Figure  12.  Geometry  for  Celestial  Path  Types  3  and  4. 

For  path  type  3,  the  apparent  zenith  angle,  8',  is  known  and  the  geometric  zenith 
angle,  8,  is  the  desired  result.  As  pointed  out  in  Chapter  2,  they  are  related 
by  the  total  bending,  4» ,  from  HI  to  the  "top"  of  the  atmosphere 

8  =  6'  +  i|>  (52) 

The  procedure  for  type  3  paths,  then,  is  to  trace  from  HI  to  the  "top"  of  the 
atmosphere  with  ANGLE  =  8'  and  calculate  8  by  Equation  (52). 

For  path  type  4,  8  is  known  at  8'  as  the  desired  result.  Equation  (52)  can¬ 
not  be  used  since  4,  is  unknown.  The  approach,  therefore,  is  to  iterate  in  a 
manner  similar  to  that  above. 

The  Newton-Raphson  Equation  is 


8 ' j+1  =  9’j  +  <e  -  8j)/T 


where  the  subscripts  refer  to  the  jth  iteration.  \e  derivative  could  be  found 
from  Equation  (52),  but  a  more  convenient  form  is  given  as  follows.  Rewrite  Equa¬ 
tion  (15)  in  terms  of  the  zenith  angles  in  Figure  12. 


Solving  for  4«  and  substituting  in  Equation  (52)  gives 


0  =  6"  +  p' 


(55) 


Differentiating 


36  38^  3_§_|_ 

36 '  36  '  36  * 


(56) 


The  second  term  was  developed  in  Section  5.2.1.  Recall  that 


(57) 


The  first  term  in  Equation  (56)  can  be  derived  in  a  manner  analogous  to  Equation 
(35)  as 


36"  _  tan  6" 
36 1  tan  6* 


6'  t  0 


(58) 


For  an  initially  horizontal  path  (S'  =  0)  we  recall  that  n  =  1.0  at  the  "top"  of 
the  atmosphere,  giving,  from  an  intermediate  stage  of  the  derivation 


36"  _  noro 

^  "  rtop  cos  6" 


where  6 '  =  0 


(59) 


This  gives 


36  36^  3_j)_!_ 

36 ’  36 '  3aQ 


(60) 


where  36"/3r'  is  given  by  Equations  (58)  or  (59)  and  3 p'/3an  is  given  in  Section 


Chapter  6 


ATMOSPHERIC  DATA  PROCESSING 


6 . 1  Atmospheric  Data  Sources 

As  outlined  in  Chapter  2,  final  processing  of  the  atmospheric  (upper-air) 
data  is  accomplished  by  the  program  BLDATM.  Output  data  form  BLDATM  may  be  input 
to  either  RAYTRA  or  the  latest  AFGL  version  of  LOWTRAN.  As  mentioned  earlier, 
the  input  upper-air  data  (which  includes  header  information)  for  BLDATM  may  come 
from  four  sources . 

a.  USAFETAC  DATSAV  raob  data  after  it  has  been  processed  first  by  ENAPRECON 
and  then  by  ENAEXTR. 

b.  Upper-air  point  analysis  data  after  it  has  been  prepared  by  the  program 
ENS2AMOD  and  processed  by  PASELECT. 

c.  Modeled  upper-air  data  derived  from  AFGL's  LOWTRAN  program  ( LOWTRN ) . 
There  are  six  models  to  choose  from  which  are  already  constructed  as  disk  files. 
They  are 


1976S . LOW; 1  -  1976  US  Standard  Atmosphere 
ARTSU . low ; l  -  subarctic  Summer  (60°N,  July) 

ARTWI . LOW; 1  -  Subarctic  Winter  (60°N,  January) 

MIDSU . LOW; l  -  Midlatitude  Summer  (45°N,  July) 

MIDWI . LOW; 1  -  Midlatitude  Winter  (45°N,  January) 

TROP I . LOW ; 1  -  Tropical  (15°N) 

The  above  models  are  profiles  that  contain  columns  of  height  (km),  temperature 
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(K),  pressure  (mb),  absolute  humidity  (g/m  ),  ozone  density  (g/m  ),  and  level 
number.  (NOTE:  Ozone  density  is  used  only  in  LOWTRN,  and  all  data  extend  upward 
to  100  kilometers  even  though  RAYTRA  does  not  use  the  data  above  50  kilometers). 

d.  Upper-air  data  which  are  input  into  BLDATM  from  the  terminal  at  time  of 
execution  (i.e.,  an  interactive  response-to-query  mode).  For  this  option,  the 
program  requires,  by  level,  input  of  height  (m,  km,  ft,  or  thousands  of  feet)  or 
pressure  (mb  or  in.  Hg)  but  not  both.  (NOTE:  A  surface  pressure  must  be  input 
when  heights  are  used,  and  a  surface  elevation  must  be  input  when  pressures  are 
used.)  Further,  it  requires  input,  by  level,  of  temperature  (°C,  °K,  or  °F)  and 
one  moisture  term.  The  moisture  terms  permitted  are  relative  humidity  (%), 
absolute  humidity  (g/m3),  mixing  ratio  (nondimensional ;  i.e.,  not  parts  per  thou¬ 
sand),  dew-point  temperature  (°C,  °K,  or  °F),  or  dew-point  depression  (°C  or  °F). 

6 . 2  Allowance  for  Missing  Data 


Missing  data  allowances  for  the  above  sources  are 


a.  Raob  data.  No  two  or  more  successive  missing  levels  or  dew-point  depres¬ 
sion  from  the  surface  to  the  level  where  the  temperature  drops  below  -40°C,  here¬ 
after  called  the  -40°C  level,  are  allowed.  (NOTE:  This  level  rarely  occurs  at 
exactly  -40°C;  therefore,  the  highest  level  falling  between  -40°C  and  -30°C  is 
used.)  All  dew-point  depressions  above  the  -40°C  level  are  considered  missing 
and  are  computed  to  the  top  of  the  sounding.  A  constant  mixing  ratio  of  3xl0-t> 
is  assumed  at  and  above  the  150-mb  level,  and  linear  height  interpolation  is  used 
to  compute  any  required  dew-point  depressions  between  the  -40°C  level  and  the 
150-mb  level.  If  a  150-mb  level  is  not  reported,  the  height  of  a  fictitious 
150-mb  level  is  computed  hydrostatically  (using  a  density  of  2.4152  kg/m  and  a 
gravitational  acceleration  of  9.7647  m/sec2,  both  taken  from  the  1976  U.S.  stand¬ 
ard  Atmosphere  [9]  for  the  150.19-mb  pressure  level).  Dew-point  depressions  for 
missing  levels  above  the  -40°C  level  are  calculated  using  Equation  (70)  below. 
All  other  data  must  be  present  and  are  assumed  correct. 

b.  Point  Analysis  data.  No  missing  temperatures,  pressures,  heights,  or 
absolute  humidities  are  allowed  from  surface  to  100  thousand  feet.  (NOTE:  BLDATM 
will  not  process  any  data  above  100  thousand  feet. ) 

c.  Modeled  data.  No  missing  data  are  allowed. 

d.  Response-to-query  data.  No  two  or  more  successive  levels  of  any  moisture 
term  are  allowed  from  the  surface  to  the  -40°C  level.  All  moisture  terms  above 
this  level  are  considered  missing.  Dew-point  depressions  are  computed  to  the  top 
of  the  sounding  by  the  same  method  described  for  raob  input.  Formulas  used  to 
convert  input  moisture  terms  into  dew-point  depressions  are  described  later  in 
this  chapter.  All  other  data  must  be  present  and  are  assumed  correct. 


6 . 3  Data  Adjustment 

All  linear  height  interpolations  use  the  following  formula 


X2  =  X1  + 


<X3  "  X1><Z2  "  V 

(z3  -  277 


(61) 


where  X2  is  the  missing  variable  at  height  Z2,  X3  and  X^  are  the  known  variables 
at  the  respective  heights  Z3  and  Z^  which  bound  Z2> 


Whatever  the  type  of  input  data,  BLDATM  will  produce  output  (for  each  obser¬ 
vation)  that  contains  header  information  an/  columns  of  MSL  heights  (km),  temper- 
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atures  (K),  pressures  (mb),  absolute  humidities  (g/m  ),  and  level  numbers.  No 
missing  data  are  allowed  in  this  output. 


Should  the  top  of  the  upper-air  data  terminate  at  a  relatively  low  level, 
such  as  850  mb,  the  data  will  still  be  processed  for  input  into  RAYTRA  or  LOWTRN. 
In  all  cases,  RAYTRA  will  perform  an  exponentially  decaying  interpolation  of 


V 


34 


refractive  moduli  from  the  top  of  the  data  to  the  50-km  level.  The  refractive 
modulus  N  at  the  50-km  level  is  assumed  to  be  zero. 

If  either  HI  or  H2  are  found  to  be  below  the  surface  level  of  atmospheric 
data,  they  will  automatically  be  adjusted  upward  to  coincide  with  the  surface 
level  by  the  program  RAYTRA. 

6.4  Pressure/Heiaht  and  Humidity  Formulas 


The  basic  formulas  used  in  BLDATM  to  compute/convert  moisture  terms  are  in 
most  meteorological  textbooks  that  deal  with  thermodynamics .  Given  the  numerous 
assumptions  in  RAYTRA  and  LOWTRN,  as  well  as  the  inaccuracies  associated  with 
most  upper-air  measurements,  these  equations  are  of  sufficient  accuracy.  They 


e  =  6.11  x  10  exp  fa)  (62) 

(Tetan's  formula  [14])  where  e  is  vapor  pressure  (mb),  Td  is  dew-point  tempera¬ 
ture  ( °C ) ,  a  is  7.5  (over  water,  T  >_  -40°C)  and  9.5  (over  ice,  T  £  -40°C),  and  b 
is  237.3  (over  water)  and  265.5  (over  ice).  The  saturation  vapor  pressure  eg  is 
computed  by  the  same  formula  except  that  the  temperature  T  (°C)  is  used  instead 
of  Td. 

w=0^2i92_e  (63) 

where  w  is  the  nondimens ional  mixing  ratio,  e  is  the  vapor  pressure  (mb)  and  P  is 
the  total  atmospheric  pressure  (mb). 

RH  =  ,64) 

s 

where  RH  is  relative  humidity  (%),  e  is  vapor  pressure  (mb),  and  eg  is  saturation 
vapor  presure  (mb) 

„„  216.494  e  lrr  , 

AH  =  - fji -  (65) 

3 

where  AH  is  absolute  humidity  (g/m  ),  e  is  vapor  pressure  (mb),  T  is  temperature 
(K),  and  216.494  is  a  constant  derived  from  the  ratio  of  the  molecular  weight  of 
water  (18.0160)  over  the  universal  gas  constant  (8.31436  x  107  erg/Mol  °K),  or 

AH  =  pW  (66) 
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where  AH  is  absolute  humidity  (g/m  ),  p  is  dry  air  density  (g/m  ),  and  w  is  the 


nondimensional  mixing  ratio. 


where  T*  is  virtual  temperature  (K);  T  is  temperature  (K),  and  W  is  the  nondimen- 
sional  mixing  ratio. 


P 


2 


exp 


(21  "  Z2] 
(T*1  +  T*2)} 


(68) 


where  P2  is  pressure  (mb)  at  height  Z2  (km),  PI  is  pressure  (mb)  at  height  Z^, 

T*,  and  T*.  are  corresponding  virtual  temperatures  (K),  g  is  gravitational  accel- 
^  J*  2 

eration  (km/sec  ),  and  R  is  the  dry  gas  constant.  (NOTE:  2g/R  is  computed  to  be 
68.283  °K  per  km  when  using  consistent  units.)  This  is  an  approximation  for  non¬ 
zero  temperature  lapse  rates.  The  error  is  negligible  for  typical  radiosonde 
data,  but  may  be  significant  for  very  thick  layers  (e.g.,  greater  than  10  km). 
It  may  be  solved  for  height,  giving 


Z2  = 


§g  <T*1  +  TV  ln  <P7> 


(69) 


Equations  (62)  and  (65)  are  combined  to  compute  dew-point  depressions  by  the 
following  equation 


DEP 


_  b  log  X 
a-log  X 


(70) 


where 


y  _  AH  (T  +  273.2) 

A  "  (216 .494 ) ( 6 . 11 ) 

and  where  DEP  is  dew-point  depression  (°C),  T  is  temperature  (°C),  and  the  other 
symbols  and  constants  are  the  same  as  in  Equations  (62)  and  (65). 

Equations  (61-70)  form  the  basis  for  all  data  computations/conversions  in 
BLDATM,  remembering  that  height,  pressure,  temperature,  and  absolute  humidity  are 
always  the  desired  output.  The  exact  form  and  sequence  of  the  equations  used 
depend  on  which  type  of  data  is  used  for  input.  Equations  (68)  and  (69)  are  used 
to  compute  pressures  from  heights  or  heights  from  pressures  whenever  arbitrary 
response-to-query  input  is  used.  Since  T*  is  a  function  of  P,  an  improved  value 
of  P2  is  obtained  by  iterating  on  Equation  (68).  Convergence  occurs  when  the 
separate  solutions  for  P2  differ  by  no  more  than  0.1  mb. 


Section  7.3  describes  a  typical  listing  of  output  from  BLDATM. 


Chapter  7 


TT 


RAYTRACE  PACKAGE  USER'S  MANUAL 


7.1  General  Overview 


The  raytrace  package  consists  of  three  main  programs:  BLDCOM,  BLDATM,  and 
RAYTRA.  The  program  BLDCOM  produces  a  command  file  that  contains  program-control 
information,  geometry  data,  and  frequency  data.  BLDATM  creates  an  atmospheric 
data  file.  These  two  files,  the  command  file,  and  the  atmospheric  data  file,  are 
the  input  files  used  by  RAYTRA,  the  raytracing  program. 

Two  versions  of  each  program  are  available.  One  version  is  designed  for 
interactive  execution  on  the  DEC10  system  at  BBNB  in  Boston;  the  other  version  is 
designed  for  batch  execution  on  the  IBM  4341  at  USAFETAC,  Scott  AFB. 

7.1.1  Purpose  of  User's  Manual.  This  chapter  will  serve  as  a  user’s  manual  for 
both  versions  of  all  three  programs.  It  contains  the  information  required  to 
successfully  execute  the  three  programs:  BLDCOM,  BLDATM,  and  RAYTRA.  Also, 
example  input  data  and  output  results  are  provided. 

7.1.2  General  Information  on  Programs.  The  program  BLDCOM,  as  its  name  implies, 
builds  (BLD)  Command  (COM)  files.  The  name  for  this  program  on  the  IBM  is 
ENABLDCO.  The  program  BLDATM  builds  (BLD)  Atmospheric  (ATM)  data  files.  The  IBM 
name  is  ENABLDAT.  The  two  files  produced  by  this  program  are  used  as  input  to 
the  RAYTRA  program.  This  is  the  actual  raytracing  program.  On  the  IBM,  it  is 
named  ENARAYTR. 

Some  symbols  that  will  be  used  throughout  this  chapter  are  explained  below. 

a.  <CR>^  is  used  to  identify  a  carriage  return  provided  by  the  user  at  the 
time  of  execution.  Obviously,  this  pertains  to  the  interactive  versions  only. 

b.  /F  is  seen  after  some  variable  names  as  they  are  requested  by  the  pro¬ 
gram.  It  signifies  that  the  variables  must  be  entered  in  real  format  with  the 
decimal  included.  Again,  this  pertains  only  to  the  interactive  versions. 

c.  /I  specifies  that  the  variable  requested  by  the  program  must  be  entered 
in  integer  format.  This  pertains  only  to  the  interactive  versions. 

Latitudes  and  longitudes  must  be  entered  as  indicated  below. 

(1)  Positive  latitude  implies  north, 

(2)  Negative  latitude  implies  south. 
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(3)  Positive  longitude  implies  west,  and 


(4)  Negative  longitude  implies  east. 

7.2  Interactive  BLDCOM  -  General 

The  BLDCOM  program  at  BBNB  interactively  creates  a  command  file  for  RAYTRA. 
This  command  file  provides  RAYTRA  with  program-control  information,  geometry 
data,  and  frequency  data. 

7.2.1  Interactive  BLDCOM  -  Structure.  BLDCOM  is  composed  of  the  main  program, 
two  system  subroutines,  and  two  user  subroutines  (Figure  13). 

—  DATE 
—  TIME 
' —  BETANS 
—  GETCON 

Figure  13.  Interactive  BLDCOM  Structure. 

The  main  program  queries  the  user  for  input  that  will  eventually  be  used  as 
the  program-control  information,  geometry  data,  and  frequency  data  in  the  ray¬ 
tracing  program,  RAYTRA. 

The  two  system  subroutines,  data  and  time,  retrieve  the  current  data  and  time 
for  header  information. 

The  two  user  subroutines,  GETANS  and  GETCON,  retrieve  user-provided  responses 
to  yes/no  questions  and  user-provided  confirmations  to  all  input,  respectively. 
Confirmations  are  required  for  all  input.  Should  the  user  key  in  the  wrong 
information,  he  can  negate  his  response  by  entering  the  character  'N1  when  con¬ 
firmation  is  requested.  Any  other  entry,  included  a  carriage  return,  will  be 
considered  a  positive  confirmation. 

7.2.2  Interactive  BLDCOM  -  Performance.  BLDCOM  is  a  simple  program.  Presently, 
it  merely  reads  the  user's  input  and  writes  it,  with  unformatted  writes  state¬ 
ment,  to  disk.  For  this  reason,  it  requires  only  three  K  code  storage  and  about 
one  CPU  second  to  handle  one  set  of  records.  Although  its  present  purpose 
appears  superfluous  (since  it  could  easily  be  an  input  routine  of  the  RAYTRA  pro¬ 
gram),  future  uses  and  requirements  justify  its  existence  as  a  separate  entity. 

7.2.3  Interactive  BLDCOM  -  Data  Base  Requirements.  Since  all  input  is  provided 
by  the  user  at  the  time  of  execution,  no  data  base  is  required  for  BLDCOM. 

7.2.4  Interactive  BLDCOM  -  Description  of  Inputs.  Other  than  informational 
questions,  the  interactive  version  of  BLDCOM  requires  a  minimum  of  four  records 


INTERACTIVE 

BLDCOM 


38 


as  input.  These  records  are  required  input  to  RAYTRA.  The  following  is  an 
explanation  of  these  input  records. 

The  first  record  is  the  label  record.  It  is  an  ASCII  description  of  the  ray- 
trace  project.  It  is  input  only  once  at  the  start  of  the  program.  It  can  be  no 
longer  than  80  characters. 


The  second  record  is  the  program-control  record.  It  consists  of  four  varia¬ 
bles:  ITYPE,  LEN,  IPRNT,  and  IOUT.  All  are  integers  and  all  must  be  separated 
by  a  comma  when  entered.  This  record  can  be  entered  an  infinite  number  of  times, 
allowing  many  combinations  of  program-control  to  be  tested  on  the  same  geometry, 
frequency,  and/or  atmospheric  data.  The  following  is  an  explanation  of  the  four 
parameters . 


ITYPE  =  1  : 

=  2  : 
=  3  : 

=  4  : 


A  horizontal  path  (not  yet  implemented). 

A  normal  path  between  two  points  in  the  atmosphere. 
Paths  to  distant  objects.  The  angle  must  be  supplied. 
It  is  the  apparent  (observed)  angle. 

Paths  to  distant  objects.  The  angle  must  be  supplied. 
It  is  the  true  (geometric)  angle. 


LEN 


=  0  :  Normal  operation  of  the  program  which  selects  the  shorter 

path  when  applicable.  (Single  segments) 

=  1  :  The  longer  path  is  selected  when  applicable.  (Two  segments) 


IPRNT  =  0 
=  1 
=  2 
=  3 


Atmospheric  data  and  level  results  are  not  printed. 
Atmospheric  data  are  printed. 

Level  results  are  printed. 

Atmospheric  data  and  level  results  are  printed. 


IOUT 


=  0  :  No  output  files  are  created. 

=  1  :  A  hard  copy  output  file  is  created  and  the  final  results, 

as  well  as  the  information  requested  by  the  variable  IPRNT, 
are  written  to  this  file. 

=  2  :  A  binary  output  file  is  created  and  the  final  results  are 

written  to  this  file. 

=  3  :  Both  the  hard  copy  and  the  binary  output  files  are  created. 


For  almost  all  purposes,  the  suggested  entry  for  the  first  record  is  2,0,3, 1. 
This  entry  will  trace  a  normal  path  between  two  objects,  cause  the  shorter  path 
to  be  taken,  cause  a  hard  copy  output  file  to  be  opened,  and  direct  level  results 
and  atmospheric  data  to  this  output  file. 

The  third  record  is  the  geometry  record.  It  consists  of  seven  parameters: 
HI,  H2,  ANGLE,  RANGE,  BETA,  H3 ,  and  DH.  All  of  these  variables  are  real  and,  as 
in  record  #1,  they  must  be  separated  by  a  comma  when  entered.  This  record  also 
can  be  entered  an  infinite  number  of  times,  thus  allowing  several  geometry  situa- 


tions  to  be  tested  against  the  same  frequency  and  atmospheric  data.  Following  is 
an  explanation  of  these  variables. 


HI 
H  2 

ANGLE 

RANGE 

BETA 

H3 


DH 


The  initial  or  sensor  altitude  in  kilometers. 

The  final  or  the  target  altitude  in  kilometers. 

The  zenith  angle  at  HI  in  degrees. 

The  apparent  or  radar  path  length  in  kilometers. 

The  Earth-centered  angle  between  HI  and  H2  in  degrees. 
The  altitude  in  kilometers  of  the  exospheric  level. 

The  RAYTRA  default  value  is  1000.0  km.  It  is  used  for 
two-segment  tangent  raytraces . 

The  height  intervals  used  between  the  50-km  level  and 
the  exospheric  level.  The  RAYTRA  default  value  is 
50.0  km. 


Since  the  purpose  of  the  raytrace  program  is  to  compute  some  of  these  variables, 
it  is  not  necessary  that  they  all  be  entered.  The  following  is  an  explanation  of 
the  possible  entries.  See  Chapter  4  for  more  detail. 


HI  must  always  be  supplied.  If  angle  is  given,  at  least  one  of  H2,  RANGE, 
BETA  must  be  supplied.  If  more  than  one  of  them  is  given,  priority,  in 
decreasing  order,  is  H2,  RANGE,  and  BETA.  If  angle  is  not  supplied,  accep¬ 
table  combinations  are 


(1) 

H2 

-  Path  assumed  tangent  at  H2  (H3  must  be  given 
in  this  case  if  the  two  segment  path  is 
desired) 

(2) 

H2 , RANGE 

-  Iterates  to  desired  RANGE 

(3) 

H2 , BETA 

-  Iterates  to  desired  BETA 

(4) 

H2 , RANGE , BETA 

-  Same  as  H2, RANGE  -  BETA  ignored 

(5) 

RANGE 

-  Horizontal  path  —  Not  implemented 

(6) 

BETA 

-  Horizontal  path  —  Not  implemented 

(7) 

RANGE, BETA 

-  Not  implemented 

Should  a  variable  not  be  supplied  by  the  user,  its  omission  must  be  indicated  by 
entering  a  -1.0  in  the  appropriate  position. 


The  final  record  is  the  ray  E-M  frequency  record.  It  consists  of  only  one 
variable,  VI.  It  is  real,  and  as  in  the  case  of  the  third  record,  the  decimal 
must  be  provided.  VI  is  the  frequency  in  wave  numbers  (CM**-1). 


7.2.5  Interactive  BLDCOM  -  Description  of  Processing.  Since  BLDCOM  only  reads 
the  user's  input  and  writes  it  to  an  output  file,  no  processing  is  required. 


7.2.6  Interactive  BLDCOM  -  Description  of  Output.  BLDCOM  produces  a  file  that 
is  written  entirely  with  unformatted  write  statements.  Following  Section  7.2.7 


(Description  of  program  execution)  is  an  example  of  the  output  produced  by 
BLDCOM. 

7.2.7  Interactive  BLDCOM  -  Program  Execution.  Figure  14  below  is  an  example  of 
a  BLDCOM  run.  It  has  been  lettered  at  various  points  to  allow  for  the  full 
explanation  that  follows  the  figure. 

The  following  is  a  step-by-step  description  of  Figure  14. 

(a)  To  start  the  execution  of  the  program,  enter  the  program  name  after  the 
BBNB  EXEC  prompt  . 

(b)  This  is  the  program  header  line.  It  includes  the  current  date,  time, 
and  version. 

(c)  This  statement  confirms  that  the  output  file  is  opened. 

(d)  The  user  can  receive  an  explanation  of  the  input  parameters  if  he  enters 
'Y'  at  this  point.  The  explanation  given  is  the  same  explanation  given 
in  Section  7.2.4.  Note  that  all  user  input  is  prompted  by  the  program 
prompt,  ' BLD> ' . 

(e)  After  every  input  by  the  user,  a  confirmation  is  required.  In  this  case 
a  carriage  return  was  entered  as  a  positive  confirmation. 

(f)  First,  the  user  is  asked  to  enter  the  label  record. 

(g)  Notice  in  this  case  that  the  label  was  not  entered  correctly,  thus  an 
'N'  was  entered  as  a  negative  confirmation. 

(h)  When  a  negative  confirmation  is  made,  the  request  is  repeated. 

(i)  Next,  the  program-control  information  is  requested. 

(j)  Next,  the  geometry  information  is  requested. 

(k)  Next,  the  frequency  information  is  requested. 

(l)  When  all  records  have  been  entered  once,  the  user  can  enter  another  set 
of  records.  As  in  the  example,  another  set  of  records  can  be  entered  if 
the  user  replies  'Y'  to  this  question.  If  an  'N'  is  entered,  the  pro¬ 
gram  will  end. 

(m)  If  another  set  of  records  is  requested,  the  user  can  change  each  record 
previously  entered  by  replying  'Y'  to  this  question.  If  an  'N'  is 
entered,  the  previous  program-control  record  will  be  repeated. 
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Figure  14.  Example  of  BLDCOM  Execution. 
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(n)  The  same  is  true  with  the  geometry  record  as  stated  above. 


(o)  Once  again,  the  same  is  true  as  above.  But,  in  this  case  the  user  has 
opted  to  change  the  original  frequency  data. 


(p)  Notice  that  if  any  response  other  than  the  permitted  'Y'  or  'N'  is 
entered,  it  is  not  accepted. 


(q)  Question  (K)  is  repeated  for  frequency  data. 


(r)  Again,  another  set  of  records  can  be  entered. 


(s)  This  is  to  confirm  that  the  output  file  has  been  closed. 


(t)  System-provided  statement  that  indicates  the  end  of  execution. 


(u)  The  CPU  and  wall  time  is  given  by  the  system. 


(v)  Control  is  returned  to  the  executive. 

The  following  is  an  example  of  the  output  produced  by  the  example  BLDCOM  exe¬ 
cution  described  above  (Figure  15).  Notice  that  the  output  file  name  is  not  user 
changeable.  It  is  RAYTRA.COM. 
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Figure  15.  Example  of  Interactive  BLDCOM  Output. 


7 . 3  Interactive  BLDATM  -  General 


The  BLDATM  program  at  BBNB  interactively  creates  an  atmospheric  data  file  to 
be  used  by  RAYTRA.  This  file  provides  RAYTRA  with  the  atmospheric  data  used  in 
the  ray tracing. 

7.3.1  Interactive  BLDATM  -  Structure.  BLDATM  is  composed  of  the  main  program, 
two  system  subroutines,  and  28  user  subroutines  (Figure  16). 

The  main  program,  BLDATM,  has  seven  segments  that  direct  the  reading  of  vari¬ 
ous  data  input  types,  processing  this  data,  performing  various  checks  on  the 
data,  and  outputting  the  data  in  the  requested  format(s). 

The  first  segments  determines  the  output  type(s)  requested  by  the  user.  The 
user  can,  at  the  start  of  the  program,  choose  to  receive  a  hard  copy  output  tile 
that  is  human  readable.  He  may  receive  a  hard  copy  of  all  observations  as  they 
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Figure  16, 


Structure  of  Interactive  BLDATM. 


are  processed  or  of  only  the  observations  that  are  discarded  during  processing. 
A  ray trace- formatted  output  file  is  always  created. 


The  six  remaining  segments  are  all  contained  within  the  main  loop  of  the  pro¬ 
gram.  They  can  be  performed  an  infinite  number  of  times.  However,  the  first 
segment  is  executed  only  once. 


The  second  segment  determines  the  input  type  to  be  processed.  The  user  can 
select  (1)  KLOW,  (2)  point  analysis  (not  yet  implemented),  (3)  model  atmospheres, 
or  (4)  interactive  input  from  the  terminal,  since  this  segment  is  within  the  main 
loop,  any  number  of  input  types  can  be  selected. 


Once  the  input  is  determined,  the  interior  loop  of  the  program  is  entered. 
The  five  remaining  segments  are  contained  within  this  interior  loop.  When  the 

•  selected  input  file  is  exhausted,  the  interior  loop  is  exited.  Control  is  then 
returned  to  the  main  loop,  and  the  user  can  select  another  input  type  or  stop 
execution. 

The  third  segment  of  BLDATM  reads  the  requested  file  one  observation  at  a 

•  time.  An  observation  is  an  upper-air  sounding  containing  height,  pressure  (un¬ 
less  input  by  the  user,  in  which  case  it  will  contain  only  pressure  or  height), 
temperature,  and  a  moisture  parameter. 


When  an  observation  has  been  read  successfully,  the  fourth  segment  directs 
•  its  processing.  If  an  error  occurs  while  reading  the  observation,  this  segment 


V 
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is  bypassed  and  the  observation  is  discarded.  While  being  processed,  various 
data  checks  are  performed  to  insure  that  the  observation  is  useful . 

If  an  observation  is  successfully  processed,  it  enters  the  fifth  segment  that 
performs  a  quick  check  to  insure  that  all  variables  required  by  raytrace  are 
available.  If  the  observation  cannot  be  processed  successfully,  the  final  check 
is  not  used.  The  observation  is  discarded. 

The  sixth  segment  is  entered  if  the  observation  successfully  completes  the 
final  check.  At  this  point,  the  observation  is  output  as  requested.  If  the  ob¬ 
servation  has  missing  data,  it  is  not  output  at  this  time.  Instead,  it  is  dis¬ 
carded. 

The  final  segment  handles  all  the  discarded  observations.  A  discarded  obser¬ 
vation  is  written  only  to  the  hard  copy  output  file,  if  one  is  requested,  since 
it  cannot  be  used  by  raytrace.  All  discarded  observations  are  labeled  on  the 
hard  copy.  Finally,  this  segment  counts  the  discarded  observations  and  keeps 
track  of  the  number  of  the  discarded  obsvations  so  the  user  can  be  given  an 
account  at  the  end  of  the  program. 

The  above  segment  concludes  the  interior  loop.  Once  all  observations  in  an 
input  file  have  been  processed,  the  user  must  select  another  input  file,  or  the 
main  loop  ends,  totals  are  written  to  the  screen,  and  the  program  ends. 

The  two  system  subroutines,  date  and  time,  retrieve  the  current  data  and  time 
for  header  information. 

Following  is  an  explanation  of  the  user  subroutines. 


(01) 

(02) 

(03) 

BLOCK  DATA 

AHTODP 

CHKDAT 

(04) 

COMPAH 

(05) 

COMPDP 

(06) 

(07) 

ERASER 

FILLAH 

(08) 

FILLDP 

(09) 

(10) 

FILMGT 

GETANS 

This  routine  is  used  to  initialize  constants. 

This  routine  converts  absolute  humidity  to  dew  point. 
Subroutine  CHKDAT  checks  to  see  that  all  data  required 
by  RAYTRA  is  available. 

COMPAH  computes  absolute  humidity  given  temperature 
and  dew-point  temperature. 

Absolute  humidities  above  the  -40°C  level  are  computed 
in  this  routine. 

This  routine  erases  ASCII  buffers. 

Missing  levels  of  absolute  humidity  are  filled  in  by 
this  routine. 

Missing  levels  of  dew-point  depressions  are  filled  in 
by  this  routir° . 

All  file  management  is  handled  by  this  routine. 

Answers  to  Y/N  questions  are  retrieved  in  this  ton¬ 
tine  . 
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(11) 

GETCON 

GETCON  retrieves  confirmations  to  the  user's  re¬ 
sponses  . 

(12) 

GETFIL 

Input  files  are  retrieved  by  this  routine. 

(13) 

GETVAL 

Integer  values  to  the  various  questions  are  retrieved 

in  this  routine. 

(14) 

HTTOPR 

Subroutine  HTTOPR  converts  height  to  pressure. 

(15) 

MRTODP 

MRTODP  converts  mixing  ratio  to  dew  points . . 

(16) 

OUTHRD 

Subroutine  OUTHRD  writes  the  hard  copy. 

(17) 

OUTRAY 

The  raytrace- formatted  output  file  is  written  by  this 

routine . 

(18) 

PROCES 

This  routine  is  responsible  for  directing  the  proces¬ 
sing  of  the  observations. 

(19) 

PRTOHT 

PRTOHT  converts  pressure  to  height. 

(20) 

RDKLOW 

KLOW  files  are  read  by  this  routine. 

(21) 

RDMODL 

Model  atmospheres  are  read  by  this  routine. 

(22) 

RDPTAN 

When  point  analysis  files  are  written  in  a  machine- 
readable  format,  they  will  be  read  by  this  routine. 

(23) 

RDTERT 

This  routine  reads  input  from  the  terminal. 

(24) 

RHTODP 

Subroutine  RHTODP  converts  relative  humidity  to  dew 
points . 

(25) 

SETFLG 

The  input- types,  when  terminal  input  is  selected,  are 
determined  by  this  routine. 

(26) 

SETKEY 

The  observations  in  KLOW  file  to  be  processed  are 
determined  by  this  routine. 

(27) 

SETPAR 

This  routine  initializes  parameters. 

(28) 

TRIMER 

This  routine  trims  input  lines  to  determine  their 
length. 

For  further  information  on  the  above  subroutines,  see  the  program  internal 
documentation. 

7.3.2  Interactive  BLDATM  -  Performance.  BLDATM  requires  about  25  K  storage.  On 
the  average,  2  CPU  second  is  required  to  process  one  observation  (an  atmospheric 
profile  of  height,  pressure,  temperature,  and  absolute  humidity). 

7.3.3  Interactive  BLDATM  -  Data  Base  Requirements.  BLDATM  accepts  three  input 
types:  (1)  input  from  KLOW- formatted  files  as  produced  by  the  ENA  extract  pro¬ 
gram  from  DATSAV ,  (2)  input  from  model  atmospheres  (the  1976  Standard  Atmosphere, 
the  arctic  winter  and  summer,  the  tropical  winter  and  summer,  and  the  tropical 
models),  and  input  provided  by  the  user  at  the  terminal.  Therefore,  to  execute 
the  program  using  KLOW- formatted  files  or  model  atmosphere-formatted  files,  a 
file  in  the  appropriate  format  is  required.  Figure  17  is  an  example  of  an  obser¬ 
vation  in  the  KLOW  format,  and  Figure  18  is  an  example  of  an  observation  in  the 
model  format. 

Notice  that  the  KLOW  observation  consists  of  a  header  line  and  a  finite  num¬ 
ber  of  data  lines. 
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Figure  17.  KLOW  -  Formatted  Observations. 


The  header  line  consists  of  the  following  variables:  The  block  station  num¬ 
ber  (BLKNUM),  the  latitude  (STALAT),  the  longitude  (STALON),  the  year,  the  month, 
the  date,  the  time  in  ZULU,  the  observation  number  in  the  KLOW  file,  the  station 
elevation  (STAELE),  the  station  type,  the  number  of  levels  in  the  observation, 
the  station  call  letters,  the  station  codes,  and  the  control  number.  In  BLDATM, 
this  header  line  is  read  in  the  format  (3X, 16, F7.2, F8. 2, IX, 312 , 14, 13 , 213, 14, 
IX, A4, 13, 12) . 

Following  the  header  line,  the  parameters  read  by  BLDATM  are:  the  level, 
pressure,  height,  temperature,  dew-point  depression,  wind  direction  and  speed, 
and  density  in  the  format  (7X, 13, F7. I, I8,2F7.1,F4.0,F5.1,E10.3) . 

For  further  explanation  of  the  parameters  used  in  BLDATM,  and  consequently  in 
RAYTRA,  see  the  internal  program  documentation. 

Note  that  model  atmosphere  is  the  1976  Standard  Atmosphere.  It,  like  the 
other  five  model  atmospheres,  consists  of  a  header  line  and  33  data  lines. 

The  header  line  consists  of  the  following:  Levels,  NM,  and  MODMSG.  Levels 
is  the  number  of  levels  in  the  obsvation,  NM  is  an  integer  that  defines  the 
observation  to  RAYTRA,  and  MODMSG  is  an  ASCII  label  that  describes  the  model. 
This  header  is  read  in  the  following  format:  ( 15, 15, 10A4) . 
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Figure  18.  Model  Atmosphere  -  Formatted  Observations. 


Each  data  line  consists  of:  Height,  temperature,  pressure,  absolute  humidi¬ 
ty,  and  density.  They  are  read  with  the  format:  (F7 . 3 , IX, F5 . 1, 3F11 .4 ) .  Notice 
that  the  variables  required  by  RAYTRA  are  in  this  type  of  file.  For  this  reason, 
no  processing  is  required  when  a  model  atmosphere  is  entered. 

7.3.4  interactive  BLDATM  -  Description  of  Input.  Along  with  several  information 
questions,  that  will  be  explained  in  Section  7.2.7  (example  execution),  three 
input  types  are  accepted  by  BLDATM. 

First,  KLOW  input  can  be  used.  This  file  is  explained  in  depth  following 
Figure  17,  and  further  information  can  be  found  in  the  program  internal  documen¬ 
tation. 

Next,  model  atmospheres  are  accepted.  Six  of  them  are  available:  The  1976 
Standard  Atmosphere,  arctic  summer  and  winter,  tropical  atmosphere,  and  the  mid¬ 
latitude  summer  and  winter  [9].  The  1976  Standard  Atmosphere  is  shown  in  Fig¬ 
ure  17.  All  model  atmospheres  follow  the  same  format.  Presently,  all  the  atmo¬ 
spheres  are  on  the  USAFETAC-DX  directory  at  BBNB  in  the  following  files 

(1)  1976S . low; 1  is  the  1976  Standard  Atmosphere, 

(2)  ARTSU. low; 1  is  the  arctic  summer  atmosphere, 
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(3)  ARTWI . low; 1  is  the  arctic  winter  atmosphere, 

(4)  TROPI . low; 1  is  the  tropical  atmosphere, 

(5)  MIDSU . low ; 1  is  the  midlatitude  summer  atmosphere,  and 

(6)  MIDWI . low; 1  is  the  midlatitude  winter  atmosphere. 

Finally,  the  user  can  enter  his  own  observations  at  the  terminal.  If  he  opts 
to  do  so,  the  program  asks  a  series  of  questions  to  determine  the  variables  and 
the  units  of  those  variables  that  will  be  entered.  At  least  three  variables  must 
be  entered.  First,  the  user  can  enter  height  or  pressure.  Height  can  be  entered 
in  feet,  thousands  of  feet,  meters,  or  kilometers.  Pressure  can  be  entered  in 
millibars  or  inches  of  mercury.  Second,  temperature  must  be  entered,  and  it  can 
be  entered  in  Celsius,  Kelvin,  or  Fahrenheit.  Finally,  a  moisture  parameter  must 
be  entered.  Five  are  available:  dew-point  temperature,  dew-point  depression, 
relative  humidity,  absolute  humidity,  or  mixing  ratio.  If  a  dew  point  is  en¬ 
tered,  it  must  be  in  the  same  units  as  temperature.  Relative  humidity  must  be 
entered  in  percent,  absolute  humidity  must  be  entered  in  grams  per  cubic  meter, 
and  mixing  ratio  must  be  entered  nondimensionally .  The  user  is  stepped  through 
the  process  of  choosing  the  input  he  desires.  First,  the  data  levels  are  en¬ 
tered.  Then,  the  header  information,  that  is,  station  elevation,  station  pres¬ 
sure,  latitude,  longitude,  and  date- time  are  requested.  Once  again,  the  user  is 
stepped  through  the  input  process. 

7.3.5  Interactive  BLDATM  -  Description  of  Processing.  BLDATM  has  many  proces¬ 
sing  routines  that  were  explained  in  Section  7.2.1.  The  primary  processing  done 
by  BLDATM  is  the  accurate  computation  of  absolute  humidity.  The  methods  used 
were  discussed  previously  is  this  technical  note. 

Secondary  processing  is  required  when  the  user  enters  an  observation  at  the 
terminal.  In  this  case,  either  height  or  pressure,  depending  upon  which  parame¬ 
ter  was  not  entered,  must  be  computed.  A  more  detailed  description  of  the 
methods  used  by  all  processing  routines  can  be  found  in  the  program  internal 
documentation . 

7.3.6  Interactive  BLDATM  -  Description  of  Output.  The  primary  responsibility  of 
BLDATM  is  to  produce  an  atmospheric  data  file  to  be  used  by  RAYTRA.  This  file 
consists  of  a  header  line  composed  of  station  elevation  (in  km),  station  pressure 
(in  mb),  latitude  and  longitude  (in  degrees),  the  time  in  ZULU,  the  date  of  the 
observation,  the  number  of  levels  in  the  observation,  and  a  description  of  the 
input  type.  Following  the  header  line  is  the  observation  itself,  composed  of 
height  in  kilometers,  temperature  in  Kelvin,  pressure  in  millibars,  and  absolute 
humidity  in  grams  per  cubic  meter.  All  output  to  this  RAYTRA  file  is  written 
with  unformatted  write  statements.  An  example  of  such  output  can  be  seen  in 
Figure  19. 
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Figure  19c. 

Figure  19.  Example  of  Output  from  Interactive  BLDATM. 

Also,  a  hard  copy,  human  readable,  output  file  can  be  written  if  the  user 
requests  it.  He  can  receive  a  hard  copy  file  of  all  the  observations  after  they 
are  processed,  a  hard  copy  of  only  the  observations  that  are  discarded  during 
processing,  or  he  can  choose  not  to  receive  a  hard  copy  at  all. 

If  an  observation  is  discarded  for  any  reason,  it  will  not  be  written  to  the 
raytr ace- formatted  output  file.  For  this  reason,  it  is  suggested  that  the  user 
opt  to  receive  at  least  a  human  readable  ouput  of  the  discarded  observation. 

7-3-7  Interactive  BLDATM  -  Program  Execution.  The  next  figure.  Figure  20,  is  an 
example  of  a  typical  BLDATM  execution.  It  is  lettered  at  various  points  and  is 
discussed  in  detail  following  the  figure. 

The  following  is  explanation  of  the  example  execution. 

(a)  To  begin  execution  of  the  program,  enter  the  program  name  after  the  exec 
prompt  . 

(b)  This  is  BLDATM  header  information. 

(c)  The  user  can  receive  status  reports  as  the  observation  is  processed. 
Note  here  that  the  confirmation  is  required  for  all  input.  An  'N'  de¬ 
noted  negative  confirmation,  and  any  other  character  denotes  positive 
confirmation  (a  carriage  return  was  used). 

(d)  The  user  can  receive  hard  copy  of  all  observations.  In  this  case,  it 
was  not  requested. 

(e)  If  hard  copy  of  all  is  not  requested,  a  hard  copy  of  the  discarded 
observation  can  be  received. 
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Figure  20a. 
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Figure  20b. 
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Figure  20d. 

Figure  20.  Example  Execution  of  Interactive  BLDATM. 

(f)  The  next  two  statements  verify  that  the  output  files  have  been  opened. 

The  output  files  are  program  prechosen  with  the  names  that  appear. 

(g)  The  input  type  is  chosen  next.  In  this  case,  KLOW  input  was  selected.  % 

(h)  The  input  file  name  is  variable  and  must  be  supplied  by  the  user. 

(i)  This  statement  verifies  that  the  input  file  was  found  and  opened. 

(j)  If  KLOW  input  is  selected,  the  user  must  decide  whether  all  observations 
in  the  file  are  to  be  processed  or  selected  observations  are  to  be  proc¬ 
essed.  In  this  case,  all  observations  were  not  processed. 

(k)  When  selected  observations  are  to  be  processed,  the  user  must  specify 
which  observations  are  requested.  Note  the  selection  must  be  sequen¬ 
tial  . 

(l)  In  this  case,  observations  2,8,53,54,55  were  selected  for  processing. 

(m)  Only  one  line  was  required  to  enter  the  observation  selection,  so  the 
observation  selection  is  complete.  Should  another  line  be  required, 
enter  an  'N'  here. 

(n)  These  statements  are  examples  of  status  reports. 
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(o)  When  an  observation  is  discarded,  the  reason  is  always  given.  Note 
that  the  discarded  observations  are  written  only  to  the  hard  copy  file, 
as  requested.  Also  note  that  the  good  observations  are  not  written  to 
the  hard  copy  file.  Instead,  they  are  written  to  the  raytrace  output 
file . 

(p)  When  the  input  file  is  exhausted,  the  user  is  informed. 

(q)  The  input  file  is  closed  when  it  is  exhausted. 

(r)  An  account  of  the  discarded  observation  is  given  as  well  as  the  number 
of  the  observation. 

(s)  The  user  can  select  another  input  type. 

(t)  This  time  the  user  has  selected  model  input. 

(u)  Once  again,  the  input  file  must  be  specified. 

(v)  Notice  in  this  case,  the  wrong  input  file  was  specified.  The  user 

negated  the  entry  with  negative  confirmation. 

(w)  The  input  file  was  requested  once  again. 

(x)  After  the  model  is  read,  the  user  can  select  another  input  type. 

(y)  Since  the  user  selected  input  from  the  terminal,  he  must  now  select  the 
variables  and  units  he  will  enter. 

( i )  Above,  pressure  was  chosen.  Here,  millibar  units  were  chosen. 

(aa)  Temperature  will  be  entered  in  Celsius. 

(ab)  Relative  humidity  will  be  entered. 

(ac)  Relative  humidity  must  be  in  percent. 

(ad)  This  is  an  explanation  to  the  user  as  to  how  the  observation  must  be 
entered. 

(ae)  After  the  final  data  line  is  entered,  a  control  Z  must  be  entered  as 
the  first  character  of  the  next  line. 

(af)  The  data  is  entered  as  shown. 

(ag)  The  header  information  is  entered  last. 
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(ah)  Once  again,  the  user  can  continue.  In  this  case,  he  has  no  more  input, 
and  the  orogram  ends . 

(ai)  A  total  of  the  discarded  observations  is  given. 

(aj)  The  hard  copy  and  ray trace- formatted  output  files  are  closed.  The  user 
can  immediately  see  this  output  by  using  the  exec  "copy"  command.  To 
get  a  printed  copy  of  the  output,  the  local  "FTP"  (file  transfer)  can 
be  used.  Also,  the  files  should  be  deleted  when  no  longer  needed. 

(ak)  A  system  account  of  CPU  and  wall  time  is  given  and  control  is  returned 
to  the  exec. 

The  following  is  an  example  of  the  raytrace-formatted  output  file.  It  is 
written  entirely  with  unformatted  write  statements.  More  observations,  in  this 
example,  will  be  written  to  this  file  than  to  the  hard  copy  output  file  since 
only  two  observations  of  seven  were  discarded. 

Note  that  all  observations,  except  the  two  that  were  discarded,  were  written 
to  this  raytrace-formatted  output  file.  The  two  discarded  observations  cannot  be 
used  by  raytrace;  however,  they  were  written  to  the  hard  copy  output  file  as 
requested  by  the  user. 

The  following  (Figure  21)  is  an  example  of  the  hard  copy  output. 

Note  that  the  two  discarded  observations  are  written  to  this  file.  Both  have 
missing  height  levels.  The  first  observation,  as  explained  in  the  status  re¬ 
ports,  was  discarded  for  this  reason.  The  second  observation,  however,  was  dis¬ 
carded  because  there  was  no  pressure  level  greater  than  or  equal  to  150  mb.  It 
was  discarded  before  the  height  check  was  made. 

7.4  Batch  BLDCOM  -  General. 

The  batch  BLDCOM  and  the  interactive  BLDCOM  differ  only  slightly.  Therefore, 
this  section  will  concentrated  on  the  differences.  The  primary  responsibility  of 
BLDCOM  is  still  to  create  the  command  file  for  RAYTRA. 

7.4.1  Batch  BLDCOM  -  Structure.  This  version  of  BLDCOM  consists  of  the  main 
program  and  one  system  subroutine  (Figure  22).  The  subroutine  VMLINE  writes  the 
USAFETAC-required  header  on  the  output. 

7.4.2  Batch  BLDCOM  -  Performance.  The  performance  of  the  batch  BLDCOM  is  the 
same  as  that  of  the  interactive  BLDCOM. 

7.4.3  Batch  BLDCOM  -  Data  Base  Requirements.  No  data  base  is  required  since  all 
input  is  made  at  the  card  reader. 
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Example  of  BLDATM  Hard  Copy  Output. 


Figure  22.  Structure  of  Batch  BLDCOM. 

7.4.4  Batch  BLDCOM  -  Description  of  Input.  The  input  parameters  in  this  version 
of  BLDCOM  are  the  same  as  those  in  the  interactive  version  of  BLDCOM.  For  this 
reason,  a  description  of  the  parameters  will  not  be  repeated  here.  The  differ¬ 
ence  is  that  in  this  batch  version  the  input  is  all  entered  at  the  card  reader. 

The  following  is  a  description  of  how  the  input  parameters  are  entered  at  the 
reader.  A  minimum  of  five  cards  is  required. 

Card  #1:  The  label  card. 

This  card  is  the  ASCII  label.  All  80  columns  are  available  as  a 
description  of  the  raytrace  project.  Card  #1  is  read  in  the 
format  (20A4). 

Card  #2:  The  program-control  card. 

Column  1:  I TYPE, 

Column  2:  LEN, 

Column  3:  IPRNT,  and 

Column  4:  IOUT. 

Once  again,  for  nearly  all  purposes,  2031  is  the  suggested  input 
for  this  card.  It  is  read  in  the  format  (411). 

Card  #3:  The  geometry  data  card. 

Columns  01-10:  HI, 

Columns  11-20:  H2 , 

Columns  21-30:  ANGLE, 

Columns  31-40:  RANGE, 

Columns  41-50:  BETA, 

Columns  51-60:  H3,  and 

Columns  61-70:  DH. 

This  card  is  read  in  the  format  (4F10.3). 

Use  a  -1.0  for  variables  not  used. 

Card  #4:  The  frequency  data  card. 

Columns  1-12:  VI. 

Card  #4  is  read  in  the  format  (F12.3). 

Card  #5:  The  recycle  card. 

Column  1:  To  enter  another  set  of  cards  4(2,  3,  and  4,  enter  a 

1 Y' .  To  end  program  execution,  enter  a  'N' . 

This  card  is  read  in  the  format  (Al). 


The  following  is  an  example  of  a  BLDCOM  input  sequence  with  two  sets  of 
records.  In  this  example,  the  suggested  card  #2  input  (2031)  is  used  in  both 
sets .  This  causes  a  normal  path  trace  between  two  points  in  the  atmosphere 
( I  TYPE  =  2),  the  shorter  path  to  be  used  (LEN  =0),  a  print  out  of  the  atmos¬ 
pheric  data  and  the  level  results  (IPRNT  =  3),  and  a  hard  copy  output  file  to  be 
created  (IOUT  =  1).  The  geometry  card  is  also  the  same  for  both  cases: 
HI  =  0.0  km,  H2  =  45.0  km,  and  ANGLE  =  60.0  degrees.  Note  that  the  other  param¬ 
eters  are  all  missing  and  are  set  to  -1.0  as  required.  The  frequency  control 
card  is  different  for  each  case.  First  a  raytrace  will  be  performed  for  VI  =  1.0 
CM**-1,  then  one  will  be  performed  for  VI  -  10.0  CM**-1. 

Card  Column: 

00000000011111111112222222222333333333344444444445555555555  . . . 
12345678901234567890123456780012345678901234567890123456789  . . . 

Raytrace  Label  Card  -  Test  Run. 

2031 

0.0  45.0  60.0  -1.0  (-1.0  for  rest) 

1.0 

Y 

2031 

0.0  45.0  60.0  -1.0  -1.0  -1.0  ... 

10.0 

N 

Note  that  the  label  record,  Card  #1  is  entered  only  once.  For  a  lull 
description  of  the  parameters,  see  Section  7.2.4. 

7.4.5  Batch  BLDCOM  -  Description  of  Processing.  As  with  the  interactive  BLDCOM, 
no  processing  is  necessary.  The  input  is  simply  read  and  written  in  binary  to  an 
output  tape. 

7.4.6  Batch  BLDCOM  -  Description  of  Output.  The  output  of  the  batch  version  of 
BLDCOM  is  a  binary  tape  that  will  be  used  as  input  for  RAYTRA.  Since  the  output 
is  binary,  an  example  of  it  cannot  be  given.  Suffice  it  to  say  that  is  is  used 
as  input  to  RAYTRA. 

7.4.7  Batch  BLDCOM  -  Program  Execution.  The  execution  of  the  batch  BLDCOM  is 
the  same  as  that  of  the  interactive  BLDCOM.  The  information  written  to  the 
printer  is  informational  and  an  example  of  the  execution  is  not  needed. 

7 . 5  Batch  BLDATM  -  General 


The  batch  version  of  BLDATM  is  very  similar  to  the  interactive  version.  The 
main  difference  is  that  input  is  originated  at  the  card  reader  and  is  controlled 
at  the  card  reader  rather  than  at  a  terminal.  For  this  reason,  only  the  differ¬ 
ences  will  be  discussed  in  this  section  of  the  technical  note. 
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7.5.1  Batch  BLDATM  -  Structure.  This  version  of  BLDATM  is  composed  of  the  main 
program,  one  system  subroutine,  and  21  user  subroutines.  The  documentation  for 
the  main  program  as  described  in  Section  7.2.1  applies  in  this  version  as  well 
(Figure  23).  The  system  subroutine,  VMLINE,  writes  the  USAFETAC-required  headers 
on  the  output.  The  user  subroutines  are  the  same,  but  seven  of  them  are  not  used 
(FILGMT,  GETANS,  GETCON,  GETFIL,  GETVAL ,  SETFLG,  and  TRIMER). 


—  ERASER 

—  CHKOAT 

—  SETKEY 

—  OUTHRO 
OUTRAY 

—  PROCES  — 

—  ROCHARD  — 

—  RDKLOW  -j 

—  R0M00L  - 

—  RDPTAN 


Figure  23.  Structure  of  Batch  BLDATM. 

7.5.2  Batch  BLDATM  -  Performance.  The  performance  of  this  version  of  BLDATM  is 
the  same  as  the  interactive  version. 

7.5.3  Batch  BLDATM  -  Data  Base  Requirements.  Data  base  requirements  for  this 
version  are  the  same  as  the  requirements  for  the  interactive  version.  The  for¬ 
mats  for  the  KLOW  and  the  model  files  are  the  same.  The  only  difference  is  that 
the  user  provided  input  is  entered  at  the  card  reader,  and  is  formatted,  rather 
than  being  entered  at  the  terminal.  This  point  will  be  expanded  in  the  next 
section. 

If  KLOW  or  model  input  is  requested,  the  user  must  insure  that  a  tape  with 
the  desired  input  is  available  prior  to  execution  (see  USAFETAC/DNE  for  tape  file 
numbers ) . 

7.5.4  Batch  BLDATM  -  Description  of  Input.  Input  for  BLDATM  is  originated  at 
the  card  reader.  In  addition  to  the  card  reader,  a  maximum  of  two  tapes  can  be 
used  as  input.  One  tape  with  KLOW- formatted  files  and  one  tape  with  model  atmo- 
sphere(s).  These  input  tapes  can  be  selected,  but  are  not  required,  since  the 
user  can  enter  an  observation  at  the  card  reader.  Following  is  a  discussion  of 
the  options  available  to  the  user. 


Card  #1:  This  card  is  required.  It  contains  two  variables:  Status  and 

HRDCOP.  Status  is  entered  in  column  1  as  'Y'  or  ' N *  .  it  an  'N' 
is  entered  in  column  1,  no  status  reports  will  be  given.  If  a 
'Y'  is  entered,  status  reports  will  be  provided  at  various  points 
throughout  the  program.  An  example  of  these  status  repot  Is  is 
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Figure  24. 


Structure  of  RAYTRA. 


provided  in  the  technical  note  following  Figure  24.  HRDCOP  is 
entered  in  column  2  as  'Y',  ’B',  or  'N'.  If  a  'Y'  is  entered,  a 
hard  copy  output  is  written  to  the  printer  after  each  observation 
is  processed.  If  a  ' B '  is  entered,  a  hard  copy  is  written  to  the 
printer  of  every  discarded  observation.  If  an  'tJ'  is  entered,  no 
hard  copy  is  printed.  This  card  is  entered  only  once  at  the 
start  of  the  program. 

Card  #1  is  read  with  a  format  of  (2A1). 

Card  #2:  This  card  is  required  to  identify  intype,  the  input  type  to  be 

processed.  It  is  entered  in  the  first  column  as 

1  =  Card  input, 

2  =  KLOW  input, 

3  =  Point  analysis  input,  or 

4  =  Model  atmosphere  input. 

Card  It  2  is  read  with  the  format  (II). 

The  following  cards  are  dependent  upon  the  input  type  selected. 

If  card  input  is  selected,  the  following  is  the  sequence  of  cards  that  follow 
Card  #2. 

Card  #A3:  This  card  contains  five  variables;  DTFLAG,  DUFLAG,  TUFLAG, 

MTFLAG,  and  MUFLAG. 

DTFLAG  is  entered  in  column  1  as: 

1  =  Height  to  be  entered,  or 

2  =  Pressure  to  be  entered. 

DUFLAG  is  entered  in  column  2  as: 

1  =  Height  entered  in  thousands  of  feet, 

2  =  Height  entered  in  feet, 

3  =  Height  entered  in  kilometers, 

4  =  Height  entered  in  meters. 
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5  =  Pressure  entered  in  millibars,  or 

6  =  Pressure  entered  in  inches  of  mercury. 

TUFLAG  is  entered  in  column  3  as: 

1  =  Temperature  entered  in  Celsius, 

2  =  Temperature  entered  in  Kelvin,  or 

3  =  Temperature  entered  in  Fahrenheit. 

MTFLAG  is  entered  in  column  4  as: 

1  =  Dew-point  temperature  entered, 

2  =  Dew-point  depression  entered, 

3  =  Relative  humidity  entered, 

4  =  Absolute  humidity  entered, 

5  =  Mixing  ratio  entered. 

MUFLAG  is  entered  is  column  5  as: 

1  =  Dew-point  depression  or  temperature  in  Celsius, 

2  =  Dew  points  in  Kelvin,  or 

3  =  Dew  points  in  Fahrenheit. 

Card  #A3  is  read  in  the  format  (511). 

Note  at  this  point,  three  variables  will  be  read  at  the  reader:  Either 
height  or  pressure,  temperature,  and  a  moisture  parameter.  If  dew-point  tempera¬ 
ture  or  depression  is  entered,  it  must  be  in  the  same  units  as  temperature. 

Card  #A4:  This  card  contains  9  variables:  STAELE,  STAPRS,  STALAT,  STALON, 

DAY,  MONTH,  YEAR,  TIMEZ,  and  LEVELS. 

STAELE  is  the  station  elevation  in  kilometers.  It  is  entered  in 
columns  1-10  in  real  format,  the  decimal  included. 

STAPRS  is  the  station  pressure  in  millibars.  It  is  entered  in 
columns  11-20  in  real  format,  the  decimal  included. 

STALAT  is  the  station  latitude  in  degrees.  It  is  entered  in  col¬ 
umns  21-30  in  real  format,  decimal  included.  A  positive  number 
should  be  entered  for  north  and  a  negative  number  for  south. 

STALON  is  the  station  longitude  in  degrees.  It  is  entered  in 
columns  31-40  in  real  format,  decimal  included.  A  positive  num¬ 
ber  should  be  entered  for  west  and  a  negative  number  for  east. 

DAY  is  the  day  of  the  observation,  an  integer.  It  is  entered  in 
columns  41-42 . 


MONTH  is  the  month  of  the  observation,  an  integer.  It  is  entered 
in  columns  43-44. 

YEAR  is  the  year  of  the  observation,  an  integer.  It  is  entered 
in  columns  45-46. 

TIME2  is  the  time  of  the  observation  in  ZULU.  It  is  entered  in 
columns  47-50. 

LEVELS  is  the  number  of  data  levels  that  will  be  entered  at  the 
reader,  an  integer.  It  is  entered  in  columns  51-52. 

Card  #A4  is  read  with  the  format  (4F1 . 3, 312 , 14, 12 ) . 

Card  #A5  through  card  #AN  is  the  data  cards,  one  for  every  level  to  be 
entered,  up  to  the  value  of  levels.  Each  card  will  contain  the  three  variables 
described  in  card  #A3.  All  variables  are  real,  and  they  are  read  with  the  format 
(3F10.3),  so  the  height  or  pressure  must  be  entered  in  columns  1-10;  temperature 
must  be  entered  in  columns  11-20;  and  the  moisture  must  be  entered  in  columns 
21-30. 

If  KLOW  input  (INTYPE  -  2  in  card  #2)  is  selected,  the  sequence  of  cards 
required  is 

Card  #B3 :  This  card  contains  the  number  of  observations  that  will  be  read 

from  the  KLOW  file.  If  all  of  the  observations  are  requested,  A 
-1  should  be  entered.  If  a  specific  number  of  observations  is 
requested,  enter  that  number. 

If  all  observations  are  requested,  no  further  cards  are  required  here;  how¬ 
ever,  if  specific  observations  are  requested,  the  numbers  of  the  requested  obser¬ 
vations  should  be  entered  on  the  following  card(s). 

Card  #B4:  This  card  contains  the  numbers  of  the  requested  observations  to 

read  from  the  KLOW  file.  It  is  read  in  the  format  (10012)  from 
1  to  the  number  of  observations  requested  (entered  in  the  pre- 
ceeding  card).  For  example,  if  05  observations  are  requested, 
05  should  be  entered  in  columns  1-2  of  card  #B3 .  Then,  card  #B4 
would  be:  0108152345.  Notice  in  this  case,  observations  #1,8,15, 
23,45  from  the  KLOW  file  will  be  read.  Also  notice  that  they  are 
entered  sequentially.  Note  that  more  than  one  card  #B4  may  be 
required. 

If  input  from  a  model  atmosphere  (INTYPE  =  4)  is  requested,  no  other  cards 
are  necessary.  Presently,  point  analysis  input  (INTYPE  =  3)  is  not  implemented.) 
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After  these  cards  have  been  entered,  the  final  cycle  card  is  required.  It 
consists  of  one  variable,  CONTNU,  and  it  should  be  entered  in  column  1.  If  a  'Y' 
is  entered,  a  new  input  type  is  required,  and  card  #2,  as  well  as  the  appropiate 
subsequent  cards,  must  be  entered.  If  an  'N'  is  entered,  the  program  ends. 

The  following  is  an  example  of  a  batch  BLDATM  runstream.  It  is  lettered  at 
various  points  and  explained  in  depth  following  the  example. 


Card  Column: 

00000000011111111112222222222333333333344444444445555555555 

12345678901234567890123456789012345678901234567890123456789 

(a)  YY 

(b)  2 

(c)  05 

(d)  0103050754 

(e)  Y 

(f)  4 

(g)  Y 
<h)  1 

(i)  25222 

(j)  .013  1013.0  28.22  177.37  010177000029 


.013 

1013.0 

28.22 

1013.0 

292.0 

4.5 

1000.0 

291.2 

6.0 

956.0 

288.6 

6.0 

921.0 

286.6 

0.6 

850.0 
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The  following  is  an  explanation  of  the  above  runstream. 

(a)  In  this  case  the  user  has  requested  status  reports  as  the  observation  is 
being  processed,  and  he  has  requested  a  hard  copy  output  of  all  observa¬ 
tion  after  they  are  processed. 

(b)  The  user  has  selected  KLOW  input. 

(c)  Five  observations  in  the  KLOW  file  are  to  be  processed. 

(d)  The  numbers  of  the  requested  observations  are:  #1,3,5,7,54. 

(e)  The  user  will  request  another  input  type. 

(f)  The  second  input  type  is  model  input. 

(g)  Another  input  type  is  requested. 

(h)  Interactive,  or  card,  input  is  requested. 

(i)  The  user  will  enter  pressure  (DTFLAG  =  2),  and  it  will  be  entered  in 
millibars  (DUFLAG  =  5).  Temperature  will  be  entered  in  Kelvin  (TUFLAG 
=  2).  The  moisture  parameter  he  will  enter  is  dew-point  depression 
( MTFLAG  =2),  and  as  required  it  will  be  entered  in  Kelvin  (MUFLAG  =2). 

(j)  The  header  line  with  station  elevation,  pressure,  latitude,  longitude, 
day,  month,  year,  timez,  and  levels. 

(k)  The  cards  that  are  in  this  section  are  the  data. 

(l)  No  more  input  is  requested,  and  the  program  ends. 

7.5.5  Batch  BLDATM  -  Description  of  Processing.  The  processing  done  in  this 
version  of  BLDATM  is  identical  to  that  done  in  the  interactive  version. 

7.5.6  Batch  BLDATM  -  Description  of  Output.  This  version  of  BLDATM  creates  an 
output  file  of  atmospheric  data  for  raytrace.  This  file,  is  written  in  binary; 
for  this  reason,  it  cannot  be  reporduced  as  an  example.  The  hard  copy  options 
are  available  as  in  the  interactive  version,  but  this  output  is  written  directly 
to  the  line  printer,  rather  than  to  a  file  or  disk.  For  this  reason,  the  status 
reports  and  the  hard  copy  output  are  combined  on  the  print  out.  However,  the 
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formats  for  all  status  reports  and  hard  copy  are  the  same  as  the  interactive 
version. 

7.5.7  Batch  BLDATM  -  Program  Execution.  As  was  mentioned,  the  status  reports 
and  the  hard  copy  are  both  written  to  the  line  printer.  For  this  reason,  the 
example  in  Section  7.2.7  is  nearly  identical  to  the  execution  produced  by  this 
version. 

7.6  RAYTRA  -  General 


The  two  versions  of  RAYTRA,  the  one  at  BBNB  and  the  one  at  USAFETAC,  are  vir¬ 
tually  identical.  Neither  are,  in  the  true  sense,  interactive.  All  control  is 
provided  by  the  two  files  produced  by  BLDCOM  and  BLDATM,  and  no  interaction  with 
the  user  is  required.  Note,  however,  that  the  programs  BLDCOM  and  BLDATM  must  be 
executed  prior  to  the  execution  of  RAYTRA. 

7.6.1  RAYTRA  -  Structure.  The  structures  of  the  two  programs  are  identical. 
The  only  difference  in  the  two  is  the  system  subroutines  used  to  provide  header 
information.  The  interactive  version  uses  date  and  time,  and  the  batch  version 
uses  VMLINE.  Other  than  this,  the  following  structure  diagram  (Figure  24)  is 
accurate  for  both  versions. 

7.6.2  RAYTRA  -  Performance.  RAYTRA  requires  about  25  K  core  for  storage.  The 
time  of  execution  depends  upon  how  many  sets  of  record  were  entered  during  BLDCOM 
execution  and  how  many  atmospheric  observations  were  entered  during  BLDATM  exe¬ 
cution.  For  one  set  of  command  records  to  be  performed  on  one  observation,  about 
one  second  of  CPU  time  is  required. 

The  methods  used  by  raytrace  can  be  found  in  the  internal  program  documenta¬ 
tion  and  in  the  previous  sections  of  this  technical  note. 

7.6.3  RAYTRA  -  Data  Base  Requirements.  The  two  input  files  required  by  RAYTRA 
are  the  two  output  files  produced  by  BLDCOM  and  BLDATM,  and  they  have  been  dis¬ 
cussed  in  detail  in  previous  sections. 

7.6.4  RAYTRA  -  Description  of  Input.  Once  again,  the  input  files  for  RAYTRA  are 
the  command  file  produced  by  BLDCOM  and  the  atmospheric  data  file  produced  by 
BLDATM.  Both  have  been  discussed  in  depth  earlier. 

7.6.5  RAYTRA  -  Description  of  Processing.  Since  the  processing  methods  of 
RAYTRA  are  the  subject  of  this  technical  note,  they  will  not  be  mentioned  here. 

7.6.6  RAYTRA  -  Description  of  Output.  RAYTRA  output  is  determined  by  the  vari¬ 
ables  IPRNT  and  IOUT,  both  entered  during  BLDCOM  execution. 

The  user,  through  his  response  is  BLDCOM  to  the  variable  IOUT,  can  choose  to 
receive  a  hard  copy  (human-readable)  output  from  raytrace.  In  the  example  execu- 


tion  of  the  interactive  BLDCOM,  Section  7 .2.7 ,  a  hard  copy  of  raytrace  output  was 
requested.  An  example  of  this  output  follows  Section  7.5.7. 

Also,  the  user,  through  input  at  the  time  of  BLDCOM  execution,  can  choose  to 
receive  a  binary  output  file.  Such  a  file  was  not  requested  in  the  sample  execu¬ 
tion  of  BLDCOM.  It  consists  of  only  the  parameters  listed  in  the  "final  results" 
section  of  the  hard  copy  output  file  (example  following  Section  7.5.7.). 

Also,  the  user  chooses  the  output  that  will  be  written  to  the  hard  copy  file. 
The  variable  IPRNT  controls  whether  levels  results  and/or  atmospheric  data  or 
neither  will  be  written  to  the  hard  copy  file.  In  the  example  execution  of 
BLDCOM,  both  were  requested  as  can  be  seen  in  the  next  section. 

Note  here  that  in  the  interactive  version,  hard  copy  is  written  to  disk  file, 
and  in  the  batch  version,  all  output  is  written  to  the  line  printer. 

7.6.7  RAYTRA  -  Program  Execution.  The  example  execution  shown  in  Figure  25  is 
the  result  of  the  interactive  version.  The  only  difference  in  the  batch  version 
is  that  the  hard  copy  output  will  be  included  in  the  status  reports.  The  example 
execution  is  lettered  at  various  points  so  that  a  full  explanation  can  follow. 

Figure  25  is  an  example  execution  of  the  interactive  RAYTRA. 
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Figure  25.  Example  Execution  of  Interactive  RAYTRA. 


The  command  file  created  by  the  example  execution  of  interactive  BLDCOM  is 
used  in  this  execution  of  RAYTRA;  however,  the  example  execution  of  interactive 
BLDATM  created  a  larger  file  than  desired  for  this  example.  Therefore,  only  the 
model  atmosphere  1976  U.S.  standard  is  used  as  atmospheric  input. 

Figure  26  is  an  example  of  the  hard  copy  output  of  RAYTRA.  It  is  the  result 
of  the  above  example  execution  using  the  input  described. 
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7.7  Example  JCL 


The  following  JCL  (Job  Control  Language)  is  required  to  execute  the  three 
program  RAYTRACE  package  on  the  USAFETAC  4341  IBM. 

//  JOB  (Job  Name  Information) 

//  ASSGN  SYS007 ,  X'TOl^X'CO'  (BLDCOM  ouput  -  RAYTRA  input) 

(If  KLOW  input  requested  through  BLDATM,  then) 

//  ASSGN  SYS008,X'T02' ,  (KLOW) 

(If  model  input  requested  through  BLDATM,  then) 

//  ASSGN  SYS009,X'T03 1 ,  (MODEL  INPUT) 

//  ASSGN  SYS010,X'T04' ,X'C0'  (BLDATM  output  -  RAYTRA  input) 

(If  point  analysis  input,  when  implemented,  requested,  then) 

//  ASSGN  SYS011,X'T05' ,X, CO" 

(If  binary  output  from  RAYTRA  requested,  then) 

//  ASSGN  SYS012 , X'T06 1 , X' CO  1 
//  TLBL  IJSY07, ' DE181S02 , BLDCODAT ' 

(If  KLOW  input  requested,  then) 

//  TLBL  IJSYS08 

(If  model  input  requested,  then) 

//  TLBL  IJSY09 

//  TLBL  IJSYS10, ' DE181502 , BLDATDAT 1 

(If  point  analysis,  when  implemented,  requested,  then) 

//  TLBL  IJSYS11 

(If  binary  output  of  RAYTRA  requested,  then) 

//  TLBL  1JSYS12 

//  Pause,  mount  tapes  for  (job  name) 

//  EXEC  DNABLDCO 

<j_Insert  data  cards  for  BLDCOM> > 

/* 

//  EXEC  DNABLDAT 

<<Jnsert  data  cards  for  BLDATM> > 

/* 

(If  KLOW  input  requested,  then) 

//  MTC  run, SYS008 

(If  model  input  requested,  then) 

//  MTC  run, SYS009 ) 

(If  Point  analysis  input,  when  implemented,  requested,  then) 

//  MTC  run, SYS011 

(If  KLOW  input  requested,  then 

//  ASSGN  SYS008 , UA 

(If  model  input  requested,  then) 

//  ASSGN  SYS009 , UA 

(If  point  analysis,  when  implemented,  requested,  when) 

//  ASSGN  SYS011,UA 
//  EXEC  DNARAYTRA 
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/* 

//  MTC  run, SYS007 
//  MTC  run, SYS010 

(If  binary  raytrace  output  requested,  then) 
//  MTC  run, SYS012 
/* 

/& 
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Chapter  8 


APPLICATION  EXAMPLES 


8 . 1  Typical  Radar  Problem 

Find  the  elevation  angle  error  and  range  error  for  a  Midway  Island  radar  that 
is  pointing  at  a  satellite  on  1  January  1977  at  1200  GMT.  The  radar  elevation 
(HI)  is  150-m  MSL,  and  its  operating  wavelength  is  10  cm  (wave  number  0.1  cm-1). 
The  measured  (apparent)  elevation  angle  of  the  radar  is  4  degrees  (zenith  angle 
is  86  degrees).  The  satellite  altitude  (H2)  is  900-km  MSL.  The  radar's  "design 
velocity"  is  c  (speed  of  light  in  a  vacuum). 

This  is  a  straightforward  example  of  the  use  of  RAYTRA.  The  first  step  for 
an  analyst  is  to  check  the  availability  of  good  upper-air  data  for  the  location 
and  time  in  question.  Since  Midway  Island  is  a  raob  site,  a  good  sounding  was 
taken.  Once  the  analyst  has  processed  this  sounding  from  the  DATSAV  data  base 
(using  ENAPRECON,  ENAEXTR  and  BLDATM  or  typing  it  into  BLDTAM  at  the  terminal), 
he  then  prepares  the  geometry  and  electrical  information  by  executing  BLDCOM  at 
the  terminal.  Finally,  he  is  ready  to  execute  RAYTRA.  After  executing  RAYTRA 
with  a  desired  printout  of  (1)  final  results,  (2)  partial  path  results,  and  (3) 
data  used,  the  analyst  will  receive  output  as  shown  in  Figures  27  through  29. 

Figure  27  indicates  a  typically  small-range  error  of  31.07  m  for  a  total 
range  in  excess  of  3000  km.  Nearly  all  this  range  error  results  from  retarda- 
tion,  not  from  curved  path  considerations.  The  elevation  angle  error,  as  indi¬ 
cated,  is  0.23640  degrees.  Figures  28  and  29  are  self-explanatory. 

8.2  Typical  Satellite  Problem 

A  satellite  at  an  altitude  (HI)  of  722  km  views  the  Earth  at  a  nadir  angle  of 
30  degrees  (zenith  angle  is  150  degrees).  The  viewing  sensor  on  board  the  satel¬ 
lite  operates  at  wavelength  of  10  microns  (wave  number  is  1000  cm-1).  What  is 
the  ground  range  from  the  point  on  Earth  directly  beneath  the  satellite  to  the 
point  on  Earth  being  viewed  by  the  satellite?  Both  points  on  the  Earth  are  at 
mean  sea  level.  A  subarctic  winter  atmospheric  model  may  be  used. 

Results  of  case  2  are  depicted  in  Figures  30  through  32.  The  ground  range, 
as  indicated  in  Figure  30,  is  425.378  km.  Note  that  the  partial  path  results 
begin  at  the  satellite's  altitude,  not  the  ground. 


Figure  27.  Case  1  Example  Output  (1st  page) 
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Figure  31.  Case  2  Example  Output  (2nd  page). 


8.3  Atypical  Space  Problem 


An  interstellar  probe,  launched  from  the  planet  Spargan  in  the  Antares  star 
system,  is  nearing  Pluto  and  beaming  "hello"  messages  toward  Earth.  A  specially 
designed  receiver,  operating  at  a  wavelength  of  10  m  (wave  number  is  0.001  cm” 1 ) , 
is  located  on  Midway  Island  and  picks  up  the  "hello"  message  at  1200  Q1T  on 
1  January  1977.  The  strongest  signal  is  found  when  the  receiver  antenna  is 
pointing  at  an  elevation  angle  of  6.4  degrees  (zenith  angle  is  83.6  degrees)  and 
at  a  bearing  of  255  degrees.  The  antenna  height  is  5-m  MSL.  What  is  the  eleva¬ 
tion  angle  error  associated  with  the  receiver  antenna? 

This  is  an  example  of  a  path  to  a  distant  celestial  object  (one  where  all 
adjacent  incoming  rays  may  be  treated  as  parallel).  The  results  are  indicated  in 
Figures  33  through  35.  As  shown  in  Figure  33,  the  total  error  is  0.16761  de¬ 
grees.  Had  this  been  treated  as  a  ray  from  one  point  to  another,  where  H2  was 
actually  at  the  default  H3  value  of  1000  km,  the  computed  angle  error  would  have 
been  the  same.  A  value  of  H2  less  than  50  km  would  have  resulted  in  a  different 
value.  Note  that  the  receiver  antenna  height  of  0.005  km  was  automatically  ad¬ 
justed  to  the  surface  elevation  (0.013  km)  of  the  Midway  Island  raob  site.  Be¬ 
cause  of  the  special  geometry  in  this  case,  the  total  angle  error  is  the  same  as 
the  total  bending. 
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Figure  33.  Case  3  Example  Output  (1st  page) . 
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Figure  34.  Case  3  Example  Output  (2nd  page). 
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Chapter  9 


LOWTRAN/FASCODE  MODELS 


The  programs  BLDTAM  and  BLDCOM  are  used  to  provide  input  to  AFGL's  latest 
version  of  LOWTRAN  (LOWTRN)  as  well  as  to  RAYTRA.  The  theory  and  instructions 
for  use  of  LOWTRN  are  covered  in  a  separate  AFGL  report  (see  USAFETAC/DNE  for 
latest  programs  codes  and  documentation).  The  input  files  for  LOWTRN  are 
LATMIN.ext  (from  BLDTAM)  and  LOWIN.ext  (from  BLDCOM).  The  upper-case  file  names 
are  mandatory,  and  the  lower-case  "ext"  is  user-determined. 

The  program  RAYTRA  is  used  to  provide  input  geometry  information  to  the  lat¬ 
est  version  of  AFGL's  FASCODE  program  (see  USAFETAC/DNE  for  latest  program  codes 
and  documentation).  The  input  file  for  FASCODE  is  FASIN.ext. 

As  of  the  writing  of  this  technical  note,  neither  LOWTRN  (latest  version 
LOWTRAN5 )  nor  FASCOD  (latest  model  version  FASC0D1)  were  fully  operational  at 
USAFETAC  on  the  IBM  4341. 

Refer  to  the  Figure  1  for  the  flowchart  that  depicts  how  LOWTRN  relates  to 
the  raytrace  program  sequence.  FASCOD  may  be  used  in  place  of  LOWTRN. 
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